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SUMMARY 


The  canonical  equations  of  motion  in  barycentric  synodical 
Cartesian  coordinates  and  momenta  are  integrable  by  means  of  recurrent 
power  series;  these  series  are  proved  to  be  convergent  for  initial 
conditions  anywhere  in  the  phase  space  except  in  the  two  phase  planes 
of  binary  collisions. 

The  integration  by  recurrent  power  series  is  extended  to  the 
variation  equations.  It  is  used  to  compute  the  monodromy  matrix 
associated  to  the  fundamental  period  of  a  periodic  orbit.  A  simple 
formula  is  derived,  which  relates  the  trace  of  the  monodromy  matrix 
and  the  characteristic  exponents. 

These  numerical  methods  are  applied  to  evaluate  the  characteristic 
exponents  of  Rabe's  Trojan  Orbits;  they  are  found  to  be  of  the  stable 
type  for  the  ovals,  and  of  the  unstable  type  for  the  horse-shoe  shaped 
orbit. 

When  the  periodic  orbit  is  symmetric  with  respect  to  the  axis 
of  syzygies,  four  Independent  variational  solutions  computed  only  over 
half  the  period  are  shown  to  be  sufficient  for  evaluating  the 
characteristic  exponents. 
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1.  INTRODUCTION 


Steffensen  (1956)  has  shown  how  the  equations  of  motion  for  the 
planar  Restricted  Problem  o  f  Three  Bodies  lend  themselves  easily  to 
the  integration  by  recurrent  power  series  in  the  time.  He  has  set  up 
the  algorithm  for  the  Lagrangian  equations  in  the  Jovicentrlc 
synodical  coordinate  system*  and  he  proved  that*  with  the  exception 
of  initial  conditions  right  at  a  binary  collision*  the  power  series 
are  convergent. 

These  formulae*  slightly  modified  on  a  minor  point*  were  used 
extensively  for  the  first  time  by  Rabe  (1961*  1962a*  1962b)  in  computing 
periodic  Trojan  orbits  and  long  period  ovals  at  in  the  Earth-Moon 

system.  A  similar  algorithm  for  tho  Lagrangian  equations  in  the 
barycentrlc  synodical  coordinate  system  has  been  proposed  by 
Fehlberg  (1964);  he  compared  it  with  a  Runge-Kutta-Nystrdm  procedure. 
Even  for  an  orbit  of  close  approach  to  one  of  the  primaries* 

Steffensen's  method  proved  Itself  more  accurate  and  time  saving. 

We  propose  here  to  adapt  Steffensen's  ideas  to  the  canonical 
equations  of  motion  in  the  barycentrlc  synodical  coordinate  system. 

This  implies  that  we  replace  this  fourth  order  system  by  a  system  of 
eight  differential  equations*  all  of  the  first  order.  In  the  same 
manner  as  Steffensen  did  for  the  Lagrangian  equations*  we  prove  that 
the  power  series  computed  by  recurrence  are  convergent  for  any  set  of 
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initial  conditions,  provided  it  does  not  belong  to  the  phase  planes 
(x  ■  -y,  y  “  0)  or  (x  ■  1  -  y,  y  *  0)  of  binary  collisions. 

Then  we  extend  the  method  to  the  variational  equations.  These 
are  linear  equations  whose  coefficients  are  functions  only  of  the 
coordinates  along  the  reference  orbit.  For  simplicity  of 
presentation,  we  think  of  computing  first  these  coefficients  in  power 
series  by  recurrence,  using  the  power  series  representing  the 
coordinates  along  the  orbit;  then  the  variational  equations  are 
integrated  in  turn  by  the  recurrent  power  series  method.  We  show  how 
the  computation  can  be  controlled  by  the  Jacobi  integral  to  be 
verified  by  each  variational  solution;  in  the  case  when  four  independent 
variations  are  computed  at  the  same  time,  a  drastic  check  is  provided 
by  verifying  at  each  step  how  the  matrix  of  this  fundamental  system 
remains  close  to  a  completely  canonical  matrix. 

Now  that  we  are  able  to  compute  accurately  and  efficiently  the 
matrlzant  (Danby  196A)  along  any  orbit  which  is  not  on  a  collision 
course  in  the  planar  Restricted  Problem  of  Three  Bodies,  we  do  not  need 
to  go  through  an  approximate  resolution  of  a  second  order  differential 
equation  of  the  Hill  type  (Darwin  1911,  Message  1959,  Rabe  1961)  when 
it  comes  to  computing  the  characteristic  exponents  of  a  periodic  orbit. 
Indeed  the  characteristic  roots  are  the  eigenvalues  of  the  matrizant  at 
the  end  of  the  fundamental  period,  a  matrix  which  is  called  by 
Wintner  (1946)  the  rr’cinouronh  matrix.  Several  elementary  properties  of 
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the  matrix  lead  to  a  simple  relation  between  its  trace  and  the  two 
non-trivial  characteristic  roots  of  the  periodic  orbit.  The  stability 
of  a  periodic  orbit  is  characterized  quite  simply  by  this  trace  Tr(T): 
if  0  <  Tr(T)  <  4,  the  characteristic  exponents  are  of  the  stable  type; 

Tr(r)  “0  or  Tr(T)  ■  4  give  the  two  indifferent  cases;  in  all  other 
circumstances,  the  characteristic  exponents  are  of  the  unstable  type. 

For  the  sake  of  completeness,  we  show  how  this  method  of 
computing  the  characteristic  roots  is  simplified  in  the  case  of 
a  periodic  orbit  which  is  symmetric  with  respect  to  the  axis  of 
syzygics.  There  we  need  to  compute  the  matrizant  over  only  half  the 
period,  and  the  characteristic  roots  are  derived  from  the  homomorphism 
axiom  satisfied  by  the  one-parameter  group  which  the  matrizant 
generates. 

These  numerical  methods  are  tested  on  Rabe's  Trojan  Orbits.  The 
initial  conditions  recorded  by  Rabe,  the  Jacobi  constants  and  the 
periods  have  been  converted  from  the  jovicentric  coordinate  system  and 
the  units  chosen  by  Darwin  into  the  barycentric  coordinate  systems  and 
the  canonical  uni  s  defined  by  Wintner.  Then  for  those  orbits  which 
need  it,  especially  the  horse-shoe  shaped  orbit,  the  initial  conditions 
have  been  improved.  The  characteristic  roots  are  computed.  For  the 
oval-shaped  orbits,  our  results  confirm  the  indication  which  Rabe  drew 
from  a  coarsely  approximate  solution  of  the  Hill's  equation  in  the 
variation  normal  to  the  orbit.  However,  for  the  horse-shoe  shaped  orbit. 
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whereas  Rabe  tentatively  suggested  a  variational  stability,  we  find 
variational  instability.  The  disagreement  has  been  confirmed  in  two 
ways:  on  our  side,  we  recomputed  Rabe's  horse-shoe  orbit  by  starting 
at  a  different  point  on  the  orbit,  thus  obtaining  another  matrizant, 
and  we  found  this  matrix  to  be  equivalent,  rwuulo  a  similitude,  to 
the  previously  obtained  matrizant.  On  his  side,  Schanzle  (1965) 
recomputed  the  Fourier  series  expansion  of  the  orbit  and  its  associated 
Hill's  equation;  he  found  in  its  normal  displacements  far  more 
oscillations  than  reported  by  Rabe.  In  fact,  his  analysis  shows  clearly 
that  in  this  case,  conclusions  drawn  from  second  order  approximate 
solution  of  Hill's  equation  are  just  meaningless.  Also,  Schanzle  applied 
the  method  devised  by  de  Vogelaere  (1950)  and  Brillouin  (1948)  for 
solving  numerically  Hill's  equation,  and  he  obtained  thereby  characteristic 
exponents  of  the  unstable  type. 

If 

2.  EQUATIONS  OF  MOTION 

The  canonical  units  for  mass,  length  and  time  are  adopted  as  they  are 
defined  by  A.  Wintner  (1946);  the  motion  of  the  planetoid  is  referred  to 
the  barycentric  synodical  coordinate  system.  Thus  the  planar  Restricted 
Problem  of  Three  Bodies  is  described  by  the  Hamiltonian  function 

(1)  3lf  =  ';(p^  +  p^)  -  (xp  -  yp  )  -  - - ^  ^ 

X  '^y  '^y  X  D,  p- 


where 
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(2a)  »  I  (x  +  p)^  +  1'^, 

(2b)  P2  *  I (x  +  p  -  1)^  +  y^l^. 


The  canonical  equations  of  motion 


admit  the  first  Integral 

(A)  C  -  -p(p  -  1)  -  2»; 

this  Jacobi  constant  C  is  so  chosen  that  C  >  3  at  the  equilateral 
equilibrium  configuration  whatever  the  mass  ratio  p  may  be. 

With  the  introduction  of  the  quantities 
(5)  R  -  (1  -  p)pj^,  S  «  PP2^. 

the  canonical  equations  (3)  may  be  replaced  by  a  system  of  eight  differential 


equations 


X  -  Px  +  y. 


y  -  Py  -  X, 


*  XX  +  yy  +  ux. 


P2P2  ■  XX  +  yy  -  (1  -  m)x. 


p^R  -  -3Rpj^, 


P  2^  "  “3Sp 2 • 


p  -  -(R  +  S)x  +  p  -  pR  +  (1  -  p)S, 

X  y 


Py  ■  *  s)y  -  Px 


which  lend  themselves  in  an  obvious  way  to  an  integration  by  recurrent 
power  series. 

The  formal  power  series 


E 

n  >  0 


x^(At)“, 


n  >  0 


(At)", 


E 


y  •  ’ 

n  >  0 


Po  “  /  ^  S  (At)  , 


E 


Px  “  ^  Pp^^')  • 
n  ^  0 


^  R„(At)“, 

n  >  0 


E 


>  «  q^(At)  , 

^  n  >  0 


n  >  0 


S  (At)' 
n 


are  introduced  into  the  differential  equations  and  coefficients  of  (At) 


are  collected  together  for  each  n»0,l,2,...  . 


In  this  manner,  for  each  n 


f 
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one  obtains  the  eight  relations 


(n+l)x  . ,  “  p  +  y  I 
n+1  n  n 


(n+l)y  . 1  “  q  -  x  , 

•^n+1  n  n 

X  ^  (p+l)r  .,r  -  X  j  (p+l)(x  ,x  +  y  y  )  +  (n+l)px  .  , 

0<7ln  0<Tin  P+1  n-p  n+1 

2  (p+i)«p+iVp  ■  „  2  <p+i)<ViVp  +  Vi^-p'  ‘ 

OlPl"  0^p_^n  t'  t'  i'  t' 


I 


0  <  p  <  n 


(p+1)  (3r  ^,R  +  r  R  )  -  0 

p+1  n-p  n-p  p+l 


X  (p+1) (3s  ^,S  +  s  S  -  0, 

-  ^  p+1  n-p  n-p  p+1 

(n+l)p  Ax(R  +S  )+q-uR+  (l-p)S  , 

n+1  o<p<n  P  "‘P  ""P  ^n  n  n 


(n+l)q  T  y  (R  +  S  )  -  p  . 

^n+1  _  -^p  n-p  n-p  '^n 

0<p<n^  ^  ^ 


Initial  conditions  evidently  give 


Xq  *  x(0),  y^  -  y(0),  p^  »  ‘^0  “  Py^^^’ 


the  definition  of  the  additional  unknowns  require  that 


r_  -  |(x„  +  +  y^''^ 


0  ''  0 


'o'  ’ 


*0  *  l<’‘o  +  “  ■  +  PqI'’  Pfl  ■  “'o^- 


-3 
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Once  the  coefficients  of  degree  0  in  ^  are  determined,  those  of 
first  degree  are  computed  from  the  formulas: 


<Vi  * 

Wi> 

+  UXj^, 

“o"i  ■ 

<Vi 

Vi> 

-  (1  -  p)x^, 

Po^i  ■ 

-^Pl^o- 

Vi  ■ 

Pi  *  ^  V  *  %  •  ""o  ■ 

’l  ■  •>'o"*o  *  ^o>  ■  Po- 


The  recurrence  steps  from  degree  n  (^l)  to  degree  n  +  1  are  given 
explicitly  by  means  of  tfie  formulas 


(6a)  (n+l)x^^,  -  +  y^, 

(6b)  (n+l)y^^,  = 

(6c)  (n+l)r„r__^j  =  2  _  <  Vl  VpVl  Vp’ Vl  Vp> 

+  (n+l)(cx^^l+x^x^^j+yj,y_^^,), 

(6d)  (n+l)s  s  =  y.  (p+l)(x  X  +y  y  -s  s  ) 

^  ^  0  <  ^n  -  1  P  ^p+Fn-p  p+1  n-p 

+  <n+l)(x„x^^j+ypy^^,-(l-i,)x_^^j), 
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(6a)  (n+l)r„R„^,  -  -  Z  (P«)  <  Vl^ 

O^p^n-1  ^ 

(6f)  (n+l)s  S  «  -  (p+l)(3s  ,S  +s  S  )-3(n+l)s  .S  , 

0  n+1  „  ^  p+1  n-p  n-p  p+1  n+1  0 

0<p<n-l 


(6g)  (n-*-l)p 


2 


X  (R  +  S  )  +  q  -  uR  +  (1  -  m)S  , 
n+1  ^  p  n-p  n-p  ^n  n  n 

0<p<n^  ^  ^ 


(6h)  (n+l)q  »  -  y  (R  +  S  )  -  p  . 

n+1  ^  -'p  n-p  n-p  n 

0<p<n^  ^ 


If  the  coefficients  at  each  step  are  conputed  in  the  order  in  which  the 
above  formulas  have  beer  written  down,  only  known  quantities  will  occur  on 
the  right  hand  side  of  the  equations. 

In  order  to  prove  that,  when  0*  power  series  are 


convergent , 

we 

introduce,  for 

every 

n  > 

1,  the  notation 

L.  . 

1 

K  * 

n 

n(n  +  1) 

and  we  show 

that,  for  every 

n  ^  2, 

the 

inequalities 

1 

^nl 

"i  n 

<  xk  r  , 

—  n 

1 

<  rk  e  , 

—  n 

1 

^nl 

"i  ^ 

IvV  • 

1 

<  sk  e  , 

~  n 

1 

Pnl 

n 

lPk„f  . 

1 

"nl 

<  Rk  e", 

~  n 

1 

^nl 

-  ‘"'n'  ’ 

1 

<  Sk 

^  n 

imply  the  inequalities 
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(7a) 


I  "i  ‘'+1 
^+11-*^+!'  • 


(7c)  k„+il  1 


(7b) 


I  ~v  ^'*■1 

'ViliyW  . 


(70)  Ivil  - 


(7f!)  iPn+ll 


(7e)  |R„+il 


(7h) 


I  -,  n+1 

"'n+lli‘'W  • 


(70  |S„^j| 


Dealing  first  with  (6a),  we  obtain 


hence  a  sufficient  condition  for  the  validity  of  (7a)  is  that 


(p  +  y)k  f"  <  (n  +  l)xk 

n  “  n+i 


Since 


n  n+2 


(n+l)k^^j^  n(n-»-l) 


1  2  2 
» "n)  ly 


for  n  >  2, 


the  more  rigid  inequality 


(8a) 


•|(p  +  y)  1  tx 


is  also  a  sufficient  condition  for  the  validity  of  (7a). 
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Treating  (6b)  in  the  same  way,  we  obtain  as  a  sufficient  condition 
for  (7b)  that 

(8b)  -|(q  +  x)  ^  ey. 

From  (6c),  we  obtain  as  a  sufficient  condition  for  (7c)  that 

(n+2)(x^+y^+?^)  2  (P+l^VlVp  -  V- 

0  ^  p  ^n-1 

But  from  the  relation 

n  +  2  =  (p  +  2)(n  -  p  +  1)  -  (p  +  1) (n  -  p) , 


we  deduce  that 


(n+2) n+1  ^p+1  ^  n-p^  n+3  ^p+2  n-p+1^  ’ 


hence  that 


n+2o 


(n+2) 


0  <  p  <  n-1 


(p+l)k  ,k  = 

tv4-1  n  — n 


p+1  n-p  n+3 


where  we  have  defined,  for  any  n  >  1, 


2  i 


o  = 

”  1  <  p  <  n  ^ 


However,  for  any  n  >  2, 


<  1  +  W  =  1  +  -2 


3  2 


3 
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so  that 


n  +  2o 

n 


f(n+l) 


and 


n  +  2o 
_ n 

n  +  3 


2 

n  +  3 


)  <  5/3. 


Consequently,  (7c)  is  verified  if 


(8c)  ^(x^  +  y^)  +  (u  +  !xqI)x  +  ly^ly  ^  (r^  -  |r)r. 

By  symmetry,  from  (6d)  we  obtain  as  a  sufficient  condition  for  (7d) 

that 

(8d)  -^(x^  +  y^)  +  ((1  -  li)  + 

We  now  address  ourselves  to  (be).  Using  the  same  relations  and  the 
same  estimates  as  in  the  preceding  case,  we  find  that 

(8e)  38^?  1  (r^  -  ■^?)R 

is  a  sufficient  condition  for  (7e),  and  by  symmetry  that 

(8f)  3S^i  1  (Sp  -  ^)S 

is  a  sufficient  condition  for  (7f). 


At  last  we  examine  (6g).  We  give  it  the  form 


(n-H)p  ,  =  -X  (R  +S  )  -  X  (R  +S  )  +  q  -  uR  +  (l-u)S  -  x  (R  +S  ) 

'^n+l  Onn  nOO  ‘n  n  n,  ,  p  n-p  n-p 

l<p<n-l'^  ^ 
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so  that  we  come  to  tlie  inequality 


(nSl)!p^  ^1  k^c  [(R  +  S)|xq|  +  (Rq  +  S^)x  +  q  +  pR  +  (1  -  u)S] 


+  x(R+S)e"  2 


1  <  p  <  n-1 


k  k 
P  n-p 


Therefrom  we  deduce  that  a  sufficient  condition  for  (7^)  is  the  inequality 


k  1(R+S)  X.  +  <R  +S  )x  +  q  +  kR  +  (l-u)S)  +  x(R+S)  ^  I,  |(  <  (n+l)k  ,pi  . 

"  0  0  0  P  ""P" 


But  the  identity 


n  +  1  »  (p+l)(n-p+l)  -  p(n-p) 


implies  that 


(n+l)k  k  -  i(i  + 

p  n-p  n  p  n-p  n+2  p+1  n-p+1 


and  hence  that 


2 


2(n-l+2o  .) 

(n+1)  ^  k  k - / 

1  **  1  p  n-p  n(n+2) 

1  p  <  n-1  ^  ' 


0  <  n  -  1, 

n-1  —  * 


However , 


14 


which  implies  that 

n-l+2o  ,  , 

- —  <  3(1  -  -)  <3. 

n  —  n 

Therefore,  a  sufficient  condition  for  (7g)  is  that 

(8g)  |((R  +  Slix^l  +  (Rq  +  S^^)x  +  q  +  pR  +  (1  -  p)S]  +  3x(R+S)  ^  pc 

By  symmetry  from  (6h),  a  sufficient  condition  for  (7h)  is  that 

(8h)  |(  (R  +  S) I y^l +  (Rq  +  S^)yl  +  3x(R+S)  £  ic . 

Now  that  we  found  sufficient  conditions  to  be  fulfilled  in  order  that 
the  inequalities  should  be  fulfilled  recurrently,  we  have  to  check  that  they 
are  compatible. 

To  begin  with,  t  can  always  be  chosen  so  large  that  (8a),  (8b),  (8g) 
and  (8h)  are  satisfied  no  matter  what  values  the  constants  possess.  Also, 
it  follows  from  (8e)  and  (8f)  that  we  must  choose 

r  <  3r^/20,  s  <  3sp/20, 

after  which  (8e)  and  (8f)  are  satisfied  provided  that  we  choose  R  and  S 
sufficiently  large.  After  this,  (8c)  and  (8d)  will  be  satisfied,  if  we 
choose  X  and  y  sufficiently  small  in  comparison  with  r  and  s.  In 
thus  choosing  small  values  for  x  and  y,  we  do  not  run  into  difficulties, 
because  the  inequalities  (7)  show  that  small  values  of  these  constants  can 
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be  compensated  by  choosing  t  sufficiently  large. 


The  inequalities  (7)  being  satisfied  for  any  n 
series  (6)  are  dominated  each  by  series  of  tlie  form 


A  +  Bt  +  C 


2 

n  >  2 


k  t 
n 


n  n 
t  , 


2,  the  recurrent 


which  is  convergent  in  the  disk  |t|  <  1/t.  Thus  the  series  (6)  are 
convergeni  in  this  disk. 


3.  VARIATION  EQUATIONS 

We  denote  by  the  vector  6  the  displacements  (6x,6y,^p  ,Ap  )  of 

X  y 

a  solution  t  ♦  (x(t)  ,y(t) ,p^(t) ,p^(t))  of  the  canonical  equations  (3). 
This  vector  is  determined  by  the  vector  variation  equation 

(9)  *  =  V(t)^. 


V(t)  is  a  4*4  matrix  function  of  the  form 


(10) 


V(t) 


-1 
.(1) 
\  f(t) 


1 

o 

.‘•(t) 

>(t) 


1 

0 

0 

-1 


wiiere  in 


(11a)  x(t)  =  -  1 

r^(t) 


(.xCO  -*•  1  _  _ ti _  r  1  -  t  (x(0  I'-l)^ 

-{(O  J  ^2^0  [  P2^t) 
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(11b)  t<(t)  -  3  +  3  l)y(t)  ^ 

Pj(t)  pj(0  P2(t)  ■  P2(t) 


The  variational  equations  (9)  are  the  canonical  system  derived  from 
the  Hamiltonian  function 


(12)  'V  ®  *j('‘^P*^  +  i‘'P*^)  ~  ('“^x^p  -  )  -  ‘j[  j(t)i‘'X^  +  2i-<(t)6x6y  +  Y(t)<5y^l. 

X  y  y  X 


Because  the  oriy.inal  Hamiltonian  function  (1)  is  conservative,  the 
equations  (9)  verify  the  .lacobi  variational  integral 


(13) 


r 


6v 


where  the  coordinates  and  momenta  in  the  partial  derivatives  should  be  given 
their  values  at  each  time  along  the  orbit. 


It  is  our  purpose  to  show  that  the  variational  equations  can  be  integrated 
by  recurrent  power  series  together  with  the  equations  of  motion. 

For  simplicitv  we  think  of  our  task  as  two-fold.  At  each  step  of  the 
recurrence,  we  first  compute  the  coefficients  in  the  power  series  representing 
i,!"  and  j.  Then  by  means  of  the  variational  equations,  we  compute  the 
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corresponding  coefficients  in  the  power  series  representing  the  displacements. 

The  recurrent  power  series  expansion  of  a,S,Y  r'^quires  several 
auxiliary  variables.  The  selection  may  vary  widely  according  to  one's 
prejudices  and  predilections  in  programming;  it  might  also  be  influenced  by 
the  kind  of  mathematical  Information  one  would  like  to  draw  on  the  side  from 
the  integration.  We  present  here  a  list  yielding  fairly  elegant  recurrence 
formulas.  Our  own  computer  program  actually  uses  a  list  with  four  less 
auxiliary  variables  needed.  Here  we  introduce  in  a  first  block 


A  » 


X  u 


B 


X 


and  in  a  second  block 


E  «  1  -  3 

G  *  AC, 

J  •  1  -  3 

F  -  1  -  3 

H  -  BD, 

K  «  1  -  3 

so  that  the  time  dependent  coefficients  in  the  matrix  V  take  the  simple 
f  orm 


Cl  -  -(RE  +  SF), 
P  *  3(RG  +  SH), 
Y  *  -(RJ  +  SK). 


We  denote 


the 


coefficients  in 


n  ■  0 


their  power  scries  in 

A  (At)*^,  and  so  on. 
n 


the  natural  way: 
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As  for  the  variations  themselves,  we  have  put 


6x 


Cn(At)". 


6v 


2 

n  >  0 


2 

n  >  0 


2 

n  >  0 


4>  (At)*^, 
n 

^j^(At)". 


For  n  ■  0,  the  coefficients  in  the  auxiliary  series  A  to  K  and  in  the 
functions,  a,6,Y  are  computed  from  the  initial  conditions  on  the  orbit,  while 
they  are  determined  in  the  variations  from  the  chosen  initial  values  for 
the  displacement.  Once  the  coefficients  of  degree  n  have  been  computed 
for  the  coordinates  and  momenta  along  the  orbit,  and  the  displacements  to 
the  orbit,  the  coefficients  of  degree  n  1  are  computed  by  meins  of  the 
formulae  (6)  to  be  followed  by  four  new  sets  of  formtilae.  The  set  that 
ought  to  be  handled  first  is 


(14a) 


r  A  ■  X  -  _ 

On  n  ^  ,  n-p  p 

0  <  p  <  :i-l  ^  ^ 


2  ^ 


A_. 


(14b) 


s  B  ■  X  -  X.  s  B  , 
On  n  ^  ,  n-p  p 

0  <  p  <  n-1 


2 


(14c) 


*  y  - 
On  n 


2  r  C. 

n-p  p 

n  ^  n  —  1  *  r 


0  <  p  <  n-1 


(14d) 


s_D 
0  n 


2 

0  ^n-1 


s  D  ; 
n-p  p 


then,  the  recurrence  should  go  through  the  formulae 
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(15a)  E  =  -3  y.  h  h 

n  „  P  n-p 

(j  ^  p  •  n 


2 


(15c)  C  =  ^  AC  , 

n  p  n-p 

0  <  p  '  n 


2 


(15c)  .1  =  -3  C  C 

n  „  ^  p  n-p 

0  ^  p  n 


(15b)  =  -3  ^  BpVp- 

O  <-  p  <  n 


2 


(15d)  H 


2 


B  D 

n  ^  P  n-p 

0  <  p  <  n 


(15f)  .  -3  ^  D^Vp- 

0  <  p  •  n 


2 


before  the  coefficients  of  the  time  functions  in  the  matrix  V  could  be 
computed  by 


(16a)  31 

n 


2 

0  ^ p  ^n 


(RE  +  S  F 
p  n-p  p  n-p 


(16b)  b 

n 


3  2 

0  ^  p  n 


(R  G  +  S  H  ), 
p  n-p  p  n-p 


(16c)  > 

n 


(R  J  +  S  K  )  . 
p  n-p  p  n-p 


Final Iv  the  coefficients  of  degree  n  in  the  variational  solutions  are 
given  bv 


(17a)  (n  +  l)f,  =  n  +  $  , 
n+1  n  n 


n-H 


(17b)  (n  +  Dm 


n 
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(17c)  (n  + 


(  »  ■'  +  r  n  )  , 

p  n-p  p  n-p 


(17d)  (n  +  l)w 

n+1 


0  <  p  <  n 


C 

p  n-p 


y  r  )  . 
P  n-p 


A  proof  of  the  convergence  for  the  series  in  the  variations 
the  same  lines  as  the  proof  we  gave  in  the  first  section  of  this 


fol lows 
paper . 


4.  CHARACTERISTIC  EXPONENTS 

We  consider  the  four  solutions  6  ^  ^  ^ ^  ^  ^  ,  6  of  the  variational 
equations  which  are  determined  respectively  by  the  initial  conditions 


x^O)  =  1, 

"y^o)  =  0, 

6p^(0)  =  0, 

^p^0'  =  0, 

y 

x^^O)  =  0. 

^  y ^  ^ (0)  =  1 , 

^p^^O)  *  0, 

'p^^O)  =  0, 

y 

x^^^O)  =  0. 

'v^^^O)  =  0. 

•p^^^O)  =  1. 

X 

>^”(0)  =  0. 
y 

IV 

X  (0)  =  0, 

IV 

‘v“(0)  =  0, 

IV 

(0)  =  0, 

IV 

‘p  (0)  =  1. 

v 

We  call  R(t;0)  the  matrix  whose  columns  are  made  of  these  four  solutions; 
this  is  nothing  else  than  the  matrizant  of  the  variational  equations  such 
t  ha  t 

K(0;())  =  I^ 

(1,  denotes  the  4  ■  4  unit  matrix.) 
k 
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The  fact  that,  for  any  t,  R(t)  is  a  symplectic  matrix,  provides 
another  check  on  the  accuracy  of  the  numerical  integration  which  produced 
the  four  fundamental  solutions.  Indeed,  the  matrix  identity 

R(t;0)J(R(t;0))'^  =  J, 

where 


is  equivalent  to  the  6  independent  bilinear  identities 

(19a)  6x^t)6y”^(t)  -  6x^^^ (t) 6y^ (t)  +  6x^ ^ ( t ) 6y t )  -  6x^''(t)6y^^(t)  «  0, 

(19b)  6x^(t)6p^^^(t)  -  6x^^^(t)6pJ(t)  +  6x^^(t)6p^''(t)  -  6x^''(t)6p^^(t)  -  1, 

(19c)  6x^(t)6p”^(t)  -  5x”^(t)6pJ(t)  +  6x^^(t)6pJ'^(t)  -  6x^''(t)6pJ^(t)  =  0, 

(19d)  6y^(t)6p”^(t)  -  6y^^^t)6p^(t)  +  fy^^t) 6pJ''(t)  -  6y^'' ( t ) 6p^^ ( t )  =  0, 

(19e)  6y^t)6p”^t)  -  6y^”  (t)  6pJ(t)  +  6y ^ ^t )  6p ( t )  -  6y^'^(t)6py^t)  «  1, 

(19f)  fip^(t)^p”\t)  -  6py\t)6p\t)  +  6py(t)6p^^/(t)  -  6p^'^(t)6py(t)  *  0. 

xy  xy  x^  xy 

It  should  be  emphasized  that  the  check  through  these  identities  is  more 
stringent  than  the  one  provided  by  the  variational  Jacobi  integrals.  Moreover, 
it  dispenses  with  checking  that  det  (R(t;0))  =  1.  For  in  a  phase  space  whose 
dimension  is  even,  the  fact  that  the  matrix  R(t;0)  is  completely  canonical 
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implies  that  its  determinant  is  equal  to  unity.  Besides,  as  it  has  been 
shown  by  Bennett  (1965),  the  matrix  R(t;0)  is  usually  ill-conditioned 
so  that  the  computation  of  its  determinant  may  become  meaningless, 
although  it  still  makes  sense  to  compute  linear  combinations  of  its  minors 
of  order  2. 

The  characteristic  roots  associated  with  the  fundamental  period  T 
of  a  periodic  orbit  are  the  roots  of  the  polynomial  equation 


(20) 


det  (R(T;0)  -  sl^)  =  0. 


It  is  known  that  it  has  at  least  two  roots  equal  to  +1,  and  that  the 
other  two  roots  Sj^  and  S2  are  such  that 


We  put 


with  j  =  ±1;  is  called  the  oharaateristio  exponent  associated  to  the 

period  T. 

But  in  the  characteristic  equation,  the  sum  of  the  roots  is  equal  to 
the  trace,  denoted  here  Tr(T),  of  the  matrix  R(T;0);  thus 

Tr(T)  =  6x^(T)  +  6y^^(T)  +  dp^^^d)  +  ap^'^d) 

X  y 

and 


(21) 


2  +  2  j  coshfi  T  =  Tr(T) . 
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4 


This  is  a  simple  expression  relating  the  characteristic  exponent  and  the 
trace  of  the  monodromy  matrix. 

Table  I  summarizes  the  discussion  concerning  the  variational  stability 
that  ensues  from  the  fundamental  relation  (21). 


TABLE  I.  VARIATIONAL  STABILITY  OF  A  PERIODIC  ORBIT  AS  DEFINED 
BY  THE  TRACE  OF  ITS  MONODROMY  MATRIX 


9. 

j 

Tr(T)  >  A 

real 

+1 

even  Instability 

Tr(T)  =  4 

0 

+1 

indifferent  case 

2  ^  Tr(T)  <  4 

purely  imaginary 

+1 

even  stability 

0  <  Tr(T)  1  2 

purely  imaginary 

-1 

odd  stability 

Tr(T)  =  0 

0 

-1 

indifferent  case 

Tr(T)  <  0 


real 


odd  instability 


Ik 


5.  CHARACTERISTIC  EXPONENTS  OF  SYMMETRIC  ORBITS 

Along  a  periodic  orbit,  the  matrix  V(t)  in  the  right-hand  member 
of  the  variation  equations  is  a  periodic  function.  Thus  the  identity 

R(t;0)  -  V(t)R(t;0) 
implies  the  identity 

R(t  +  T;0)  -=  V(t)R(t  +  T;0). 

Accordingly,  since  R(t;0)  is  the  matrlzant  of  the  variation  equations, 
we  have  the  identity 

R(t  +  T;0)  *  R(t;0)R(T;0); 
in  particular,  at  t  *  -T/2, 

R(T/2;0)  -  R(-T/2;0)R(T;0). 

Therefore 

R(-T/2;0)[R(T;0)  -  sl^]  -  R(T/2;0)  -  sR(-T/2;0), 

and  this  proves  that  the  characteristic  equation  (20)  is  equivalent  to  the 
equation 

(22)  det[R(T/2;0)  -  sR(-T/2;0)]  -  0. 

The  equivalence  just  stated  is  valid  for  any  periodic  orbit,  whether 
symmetric  or  asymmetric. 
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We  shall  now  consider  the  particular  case  of  a  symmetric  periodic 
orbit . 

In  order  that  the  orbit  t  ->■  (x(t)  ,y  (t )  ,p^  (t)  .p^ (t)  )  be  symmetric 
with  respect  to  the  syzygy  axis  Ox,  it  is  necessary  and  sufficient  that, 
at  any  time, 

x(-t)  =  +  x(t),  y(-t)  =  -y(t). 

Consequently,  along  such  an  orbit, 

a(-t)  =  a(t), 

6(-t)  *  -6(t), 

Y(-t)  -  Y(t). 

Thus,  if  we  put 

10  0  0 
0-100 
0  0-10 
0  0  0  1 

the  fact  that  the  orbit  is  symmetric  with  respect  to  Ox  implies  in 
the  variation  equations  that 

V(-t)  -  -SV(t)S. 

Accordingly,  the  substitution  t  -►  -t  in  the  identity 

R(t;0)  =  V(t)R(t;0) 
provides  the  Identity 


R(-t;0)  =  SV(t)SR(-t;0) 
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2 

or,  since  S  *  the  identity 

^(SR(-t;0)]  -  V(t)(SR(-t;0)). 

Hence,  using  again  the  fact  that  R(t;0)  is  the  natrizant  of  the  variation 
equations,  we  come  at  last  to  the  identity 

(23)  SR(-t;0)  *  R(t;0)S. 

In  particular,  at  time  t  *  T/2, 

SR(-T/2;0)  *  R(T/2;0)S. 

Therefore,  the  determinantal  equation  (20)  takes  the  form 

(24)  det(R(T/2;0)  -  sSR(T/2;0)S]  -  0; 

it  proves  that  the  computation  of  the  characteristic  exponents  for  a 
symmetric  orbit  requests  the  integration  of  the  variational  equations  over 
only  half  a  period. 

This  proposition  has  been  stated  first  by  Moulton  (1914);  but  the 
proof  he  gives  for  it  depends  too  closely  on  the  particular  problem  he  is 
dealing  with,  namely  the  orbital  stability  of  Jupiter's  satellite  VIII,  and 
it  is  incorrect  on  several  points,  de  Vogelaere  (1950)  has  shown  how  to 
use  it  for  extracting  numerically  from  Hill's  equation  the  characteristic 
exponents  of  a  symmetric  orbit.  We  now  propose  to  do  the  same  for  the 
solution  of  equation  (24). 
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First  we  observe  that  the  matrix  identitv  (2i)  is  equivalent  to  the 
scalar  identities: 


6x^(-t)  = 

•^x^  (t)  , 

11/ 
■X  ( 

-t)  = 

-'X^^t), 

■> 

< 

1 

II 

-^y^t), 

X  11/ 
Ay  ( 

-t)  = 

w^^t), 

II 

1 

X 

c. 

-  1/  ^ 
-rtp^(t) , 

‘P”< 

-t)  = 

11/  X 

^Px  (t). 

^p^(-t)  = 

y 

'^p^  (t ) , 

y 

-t)  = 

,  11,  , 

-  f  P  (  t  )  , 

y 

•6x  ^  t ) , 


‘'X  (-t) 


X  IV,  , 
6x  (t ) , 


^y‘^-t)  - 


■'y'^V). 


*p‘‘'(-0  . 


.pf  (-t)  = 


.111/  X 

(-t)  = 


X  III/  ^ 
-^P..  (O, 


IV,  - 
P  (-t)  = 
V 


.  IV/  V 

P  ( t ) 
V 


Therefore,  the  determ inanta 1  equation  (24)  can  be  written  explicit Iv  as 


i 

1 

(l-s)^x^T/2) 

(l+s)Ax^^T/2) 

(l+s)‘x^^^T/2) 

(l-s).^x^^'(T/2) 

(l+s).^y^T/2) 

(l-s).S-^^T/2) 

(1-s)  w^^\r/2) 

(l+s)^v'^’(T/2) 

(25) 

(l+s)^p^T/2) 

X 

(l-s)^p^^T/2) 

X 

(l-s)‘p^^^T/2) 

X 

(1+s)  p'^  (  I7J) 

X 

(l-s)‘p^T/2) 

y 

(l-f-s)*p^NT/2) 

y 

(l+s)^p^^^T/2) 

y 

(1-s)  p‘^'(i72) 

We  put 

(l,2)(3,4)  =  l  ^x\T/2)^y 

IV  I  \’ 

(T/2)-^x  (T/2)^v 

^T/2)] 

rx  II 

(T/2)  ‘p^^^T/2)-‘p’  ^  ^ 

V  X 

(T/2)  pj/  (T/2)  1  , 
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(1.3) (2,4)-[6x^(T/2)6p^'^(T/2)-^x^''(T/2)6p^(T/2)] 

X  X 

(5y^^T/2)6Py^^T/2)-6y”^(T/2)(Sp“(T/2)l  , 

(1.4) (2.3)*(6x^T/2)6pJ'^(T/2)-6x^'^(T/2)6Py(T/2)) 

[6y^^T/2)6py^T/2)6y^“(T/2)6p^^T/2)], 

X  X 

(2.3) (l,4)-(6y^T/2)6p^^(T/2)-6y^'^(T/2)6pJ(T/2)l 

X  X 

[6x^^(T/2)6Py^^(T/2)-6x^”(T/2)(SpJ^(T/2)]  , 

T  TV  TV  T 

(2.4) (l,3)  =  (6y'(T/2)6py^T/2)-6y"''(T/2)6p\T/2)l 

[6x^^(T/2)6p^^^(T/2)-6x^^^(T/2)6p^^(T/2)] , 

(3.4) (l,2)*(6p^(T/2)6p^^(T/2)-6p^'^(T/2)6p^T/2)) 

XV  X  y 

[6x^^(T/2)6y^^^(T/2)-6x“^(T/2)6y^^(T/2)) 
so  that  the  determinantal  equation  (25)  takes  the  form 

(26)  (l+s)^(2,3)(l,4)+(l-s)^(l,4)(2.3)+(l-s^)^ 

[(1,2)(3,4)-(1,3)(2.4)-(2,4)(1,3)+(3,4)(1,2))=0 

In  this  form  we  see  that,  in  order  for  the  characteristic  equation  to  have  a 
root  equal  to  +1,  it  is  necessary  and  sufficient  that 
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(27)  (2.3)(1,4)  =  0. 

Moreover,  the  matrix  R(T/2;0)  has  a  determinant  equal  to  +1.  This 
determinant  being  obtained  by  making  s  ■  0  in  (26),  we  obtain  the 
relat ion 


(28)  1  -  (1.4)(2,3)  «  (1,2)(3.4)  -  (1,3)(2,4)  -  (2,4)(1,3)  +  (3,4)(1,2). 

In  view  of  (27)  and  (28),  the  characteristic  equation  (26)  takes  the 
simple  form 

(l-s^)l(l-s)^l,4)(2,3)  +  (l+s)^[l-(l,4)(2,3))f  •  0. 

Consequently,  the  non-trivial  character ist Ic  roots  are  the  solutions  of  the 
quadratic  equation 

1  +  2(l-2(l,4)(2,3))s  s^  -  0. 


As  in  the  general  case,  we  write  these  two  roots  as 


s 


1 


with  j  -  ll,  so  that 


(29)  1  +  J  coshCT  -  2(1,4)(2,3). 


There  ensues  from  it  information  concerning  the  stability  of  a  symmetric 
periodic  orbit;  the  conclusions  are  summarized  in  Table  II. 
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TABLE  II.  VARIATIONAL  STABILITY  OF  A  SYMMETRIC  PERIODIC  ORBIT 
AS  DEFINED  BY  ITS  MONODROMY  AT  HALF  THE  PERIOD 


D 

i 

(1,4)(2,3)  >  1 

real 

+1 

even  instability 

(1,4)(2,3)  -  1 

0 

+1 

indifferent  case 

1/2  1  (1,4)(2,3)  <  1 

purely  imaginary 

+1 

even  stability 

0  <  (1,4)(2,3)  1  1/2 

purely  imaginary 

-1 

odd  stability 

(1,4)(2.3)  -  0 

0 

-1 

indifferent  case 

(1.4)(2,3)  <  0 

real 

-1 

odd  instability 

6.  RABE'S  TROJAN  ORI.ITS 

In  order  to  numerically  calculate  the  orbits  described  in  this  section 
a  double  precision  FORTRAN  IV  program  for  computing  asymmetric  periodic  orbits 
in  the  plane  restricted  problem  of  three  bodies  was  run  on  the  IBM  7094  computer. 
In  addition  to  calculating  the  basic  dependent  variables  of  an  orbit  with  given 
initial  conditions,  the  program  provides  for  the  calculation  of  the  four  in¬ 
dependent  solutions  and  of  the  variational  equations.  If 

(for  a  given  value  of  the  period)  the  initial  values  are  such  that  the  orbit 
is  truly  periodic,  the  characteristic  roots  of  the  periodic  orbit  may  be 


calculated.  If  the  orbit  is  not  as  close  to  being  periodic  as  may  be  desired 


31 


(i.e.,  if  either  | x(T)  -  x(0) | ,  |y(T)  -  y(0)|,  |p^(T)  -  p^(0)[,  or 

|Py(T)  -  p  (0) I  are  larger  than  some  given  constant,  such  as  for  example 
-8  -10 

10  or  10  ),  then  the  program  provides  for  changing  the  initial 

conditions  slightly  so  that  an  orbit  will  1,  obtained  which  is  closer  to 
being  periodic  with  the  given  period  T.  Of  course,  reasonable  approxim.at  ions 
to  the  initial  conditions  are  necessary;  it  is  not  expected  that  the  program 
will  be  asked  to  find  a  periodic  orbit  of  given  period  starting  with  arbitrary 
guesses  for  initial  conditions. 

In  order  to  explain  this  method  of  "improving  orbits  to  make  them  truly 
periodic"  some  simplifying  notation  will  be  used.  Let  x(t*XQ)  represent 
the  vector  whose  four  components  are  x(t;x  ,y  ,p  ,p  ),  y(t;x  ,y  ,p  ,p  ), 

p  (t;x  ,y  ,p  ,p  '>  .nd  P  (t;x  ,y  ,p  ,p  ),  where  these  four  variables 
u  u  Xq  y  u  u  Xq  y^ 

represent  the  soluticn  of  equations  (3)  with  the  particular  initial  conditions 

-*■1  ->  ->11  ^  *111  »  -►IV  -*• 

(xQ,yQ,p^  ,P^  ).  Similarly,  X  (t;Xp),  6  (t;xQ),  i  6  represent 

the  feur  linearly  independent  solutions  of  the  variational  equations.  For  the 
given  period  T,  it  is  assumed  that 


x(T;xq)  ^  x(0;xq) 


alth.'ugh  these  vectors  are  not  "too  far"  from  being  equal.  Considering  the 
first  crdei  variations,  it  is  desired  to  consider  solu:ions  of  the  form 

MO)  x(t;xQ  ^Xy)  =  x(t;xQ)  <5^ 

^  Jll,  ^  ^  ,IV,  -«■  , 

+  “‘a*  *"’^0 


where  Ax^  has  components  Ax^,  Ay^,  ‘.p  ,  Ap 

^0  0  x^  y^ 
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Setting  t  »  0  and  remembering  the  initial  conditions  (18),  it  is  seen 
that  one  must  have 


(31) 


Ax 


0 


Ay 


0 


Ap_ 


0 

0 


If  it  is  decided  to  use  the  given  period  T  as  a  fixed  quantity,  then  it 
will  be  desired  to  have  the  solution  satisfy  the  conditions 


x(T;Xq  +  AXj^)  *  ^  ^0 


Substituting  this  and  equations  (31)  Into  equations  (30)  gives  the  results 


(32)  Xq  AXq  =  x(T;Xq)  AXq5^(T;Xq)  +  Av^6^^(T;Xq)  +  Ap^  6  ^  ^  ^  (T;Xq)  +  Ap^  6^'^(T,Xp) 

These  linear  equations  are  solved  for  the  components  of  the  unknown  %•  If 
the  new  guesses  Xq  ^Xq  initial  conditions  still  do  not  produce  a 

satisfactory  periodic  orbit,  the  process  may  be  repeated.  There  is,  however,  a 
limit  on  how  extremely  close  to  a  periodic  orbit  one  can  come  by  means  of  this 
iterative  procedure.  In  general,  we  found  that  with  16-place  arithmetic  we 


could  ensure  that  no  component  of  x(T»Xq)  would  differ  from  the  corresponding 

— 10 

component  of  x(0»Xq)  by  more  than  10  .  However,  if  much  better  initial 


guesses  are  known  and  used,  the  equations  (32)  become  extremely  ill-conditioned. 
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If  we  had  an  exact  periodic  orbit,  would  equal  x(T;x^),  and 

equations  (32)  would  become  precisely  equations  (20)  with  s  *  +1.  That 
is,  if  we  try  to  solve  equations  (32)  under  these  conditions,  it  means  that 
we  are  trying  to  find  an  eigenvector  (corresponding  to  the  known  eigenvalue 
+1)  for  the  matrix  R(T;0). 

In  describing  his  Trojan  orbits,  Rabe  uses  a  jovicentric  synodical 
coordinate  system.  The  mass  of  the  sun  is  taken  as  unity  and  the  period  of 
Jupiter  is  2^  •^1-u .  If  we  use  asterisks  to  represent  variables  used  by  Rabe, 
the  transformation  equations  giving  our  units  in  terms  of  Rabe’s  are 


(33) 


dx 

dt 

In 

dt 


1  -  u  -  X* 


-y* 


. 

dt* 

dt* 


1 


'1-u 


T* 


C  -  (l-p)C*. 

The  orbital  values  given  in  Tables  I  and  II  of  Rabe  (1961)  and  Tables  I  and  II 
of  Rabe  (1962)  were  transformed  by  means  of  these  equations  and  by  the  equations 


dx 

Px"-  y 


dy 

P  ■  +  X 

y  dt 
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to  give  the  initial  values  in  Table  III  and  the  values  of  T  and  C  given 
in  Table  IV. 

In  each  case,  the  period  T  as  given  in  Table  IV  was  taken  as  a  fixed 
parameter,  and  a  more  accurate  periodic  orbit  was  calculated  as  described  in 
the  preceding  paragraphs.  For  these  new  orbits,  the  characteristic  roots 
were  calculated.  In  Table  V  are  given  the  new  initial  values  for  the  periodic 
orbits.  Table  VI  lists  the  Jacobi  constant  as  well  as  the  trace  of  the 
monodromy  matrix  and  the  characteristic  roots  corresponding  to  the  periodic 
orbit . 
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TABLE  III. 

Rabe  ’  s 
Parameter 

STARTING  VALUES 

FOR  PERIODIC  ( 

ORBITS— 

(to  the  units  of 

transf  ormed 
this  paper. 

d 

o 

^0 

^0 

'Px>0 

‘Py'o 

1.0025 

.50029612A643 

-.868190470 

.864939351691 

.498421029587 

1.0050 

.501546124643 

-.870355530 

.863857161536 

.497798523295 

1.0075 

.502796124643 

-.872520600 

.862778849534 

.497178495820 

1.0100 

.504046124643 

-.874685658 

.861704291736 

.496561161060 

1.0125 

.505296124643 

-.876850720 

.860633590098 

.495946217159 

1.0150 

.506546124643 

-.879015780 

.859566642668 

.495333906001 

1.0175 

.507796124643 

-.881180850 

.858503443453 

.494724073660 

1.0200 

.509046124643 

-.883345912 

.857444086401 

.494116614187 

1.0225 

.510296124643 

-.885510980 

.856388739438 

.493511365658 

1.0250 

.511546124643 

-.887676040 

.855340133254 

.492908300088 

1.0275 

.512796124643 

-.889841100 

.854287389045 

.492306843750 

1.0300 

.514046124643 

-.892006166 

.853245632595 

.491707425438 

l.OAOO 

.519046124643 

-.900666420 

.849111400136 

.489392102889 

1.0500 

.524046124643 

-.909326674 

.844638938108 

.487334093530 

1.0600 

.529046124643 

-.917986928 

.840686966660 

.484695798129 
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TABI.E  IV.  PERIODS  AND  JACOBI  CONSTANTS  (using  initial  conds.  in  Table  III) 


Rabe's  Parameter 

Period 

Jacobi  Constant 

“o 

T 

1.0025 

78.11824642410334 

3.0000046197401 

1.0050 

78.2124113456173* 

3.0000184246002 

1.0075 

78.37095698004493 

3.0000413337195 

1.0100 

78.59711486899776 

3.0000732673572 

1.0125 

78.89553723182391 

3.0001141457676 

1.0150 

79.27284722811138 

3.0001638928252 

1.0175 

79.73765896723368 

3.0002224296547 

1.0200 

80.30133787108203 

3.0002896819513 

1.0225 

80.97856094134159 

3.0003655816422 

1.0250 

81.79018812928275 

3.0004502334866 

1.0275 

82.75360773001720 

3.0005428247134 

1.0300 

83.91284074382773 

3.0006443168953 

1.0400 

91.82099333962959 

3.0011363774494 

1.0500 

117.7920128526034 

3.0017245907688 

1 . 0600 

183.5145658929861 

3.0024430527304 
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TABLE  V.  STARTING  VALLES  FOR  PERIODIC  ORBITS  (PRESENT  RESULTS) 


Rabe's  Parameter 

“o 

^0 

^0 

(approximate) 

1.0025 

.  500294701009 

-.868189394186 

.864940407774 

.498420857938 

1.0050 

.501546202315 

-.870355763449 

.863857082677 

.497798432577 

1.0075 

.502795932420 

-.872520424844 

.862778961872 

.497178566762 

1.0100 

.504046133194 

-.874685581794 

.861704299953 

.496561215540 

1.0125 

.505296118194 

-.876850529922 

.860633607448 

.495946396251 

1.0150 

.506546050763 

-.879015868171 

.859566684583 

.495333815576 

1.0175 

.507796128489 

-.881180842237 

.858503483592 

.494724052434 

1.0200 

.509045980033 

-.883346025885 

.857444171146 

.494116467086 

1.0225 

.510296189709 

-.885510985492 

.856388586745 

.493511547160 

1.0250 

.511545765989 

-.887676141744 

.855340268383 

.492908273233 

1.0275 

.512795779277 

-.889841341844 

.854287586217 

.492306477874 

1.0300 

.514045824807 

-.892006380486 

.853245786222 

.491707120286 

1.0400 

.519045209040 

-.900666970996 

.849111898087 

.489391219068 

1.0500 

.524045017193 

-.909327242850 

.844639557053 

.487333137683 

1.0600 

.529053778602 

-.917980837710 

.840684004699 

.484703533306 
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TABLE  VI.  JACOBI  CONSTANTS,  TRACES,  AND  CHARACTERISTIC  ROOTS  USING 
INITIAL  CONDITIONS  FROM  TABLE  V  AND  WITH  SOME  PERIODS  AS 
IN  TABLE  IV 


Rabe's  Parar'eter 

Jacobi  Constant 

Trace 

Characteristic  Roots 

‘^O 

(approximate) 

C 

1.0025 

3.0000046137301 

.43885321 

-.78057339 

1  .625064141 

1.0050 

3.0000184264690 

.32970410 

-.83514795 

t  .550025371 

1.0075 

3.0000413312188 

.17882263 

-.91055869 

1  .413313741 

1.0100 

3.0000732665714 

.04109871 

-.97945065 

i  .201683991 

1.0125 

3.0001141435380 

.00832747 

-.99583627 

t  .091159931 

1.0150 

3.0001638941519 

.21110820 

-.89444590 

1  .447176171 

1.0175 

3.0002224306764 

.79229906 

-.60385047 

t  .797097621 

1.0200 

3.0002896834712 

1.81711725 

-.09144138 

1  .995810461 

1.0225 

3.0003655812550 

3.09222181 

.54611090 

t  .837712891 

1.0250 

3.0004502346940 

3.96383295 

.98191648 

t  .189314641 

1.0275 

3.0005428266252 

3.44565407 

.72282704 

t  .691029001 

1.0300 

3.0006443182264 

1.33281279 

-.33359361 

t  .942716981 

1 . 0400 

3.0011363784399 

.11223160 

-.94388420 

t  .330276581 

1.0500 

3.0017245910778 

1.17049332 

-.41475334 

t  .909933881 

1.0600 

3.0024430503909 

4.07526522 

1.3145467  and. 76071849 

0  IB  1/4  3/S  1/3  )/S  3/4  7/%  I 


*  I 

Fig.  1.  Variation  t'x  for  Trojan  Orbit  (C  =  3 . 0000603664A98  and 
T  *  78.50504948179581)  versus  time  (the  unit  of  time  is  the 
period  T) . 
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The  orbits  described  by  the  parameters  in  Tables  and  VI  are 

practically  the  same  as  Rabe  computed.  The  only  reason  we  refined  them 

slightly  using  our  methods  and  lb-place  arithmetic  was  to  enable  us  to 

calculate  the  characteristic  roots  with  good  accuracy;  this  point  will  be 

discussed  later.  The  fact  that  the  entries  in  Table  V  are  quite  different 

from  corresponding  entries  in  Table  III  does  not  mean  that  the  orbits 

themselves  are  very  different;  our  iteration  method  has  merely  allowed  for 

starting  at  a  slightly  different  point  on  the  orbit.  The  only  real  difference 

between  the  orbits  is  due  to  the  fact  that  we  took  Rabe's  5-decimal  place 

value  for  the  period  as  being  an  exact  constant,  and  found  the  more  accurate 

orbit  which  has  this  exact  period.  Another  indication  that  there  is  not  too 

much  difference  between  Rabe*s  original  orbits  and  our  modified  ones  is  that 

the  Jacobi  constants  as  given  in  Table  IV  do  not  differ  from  the  new  Jacobi 

-9 

constants  as  given  in  Table  VI  by  more  than  6  «  10 

In  the  way  in  which  our  computer  p' ogram  was  used  for  the  cases  represented 
in  the  tables,  the  number  of  terms  of  the  series  was  fixed  at  from  16  to  20. 

The  step  size  was  then  determined  so  that  none  of  the  last  three  terms  in  any 
of  the  series  would  be  more  than  10  At  the  end  of  the  orbit,  the  number 

of  terms  taken  in  the  series  might  in  some  cases  be  diminished  in  order  to  keep 

underflow  from  occurring.  We  considered  an  orbit  to  be  periodic  if  no 

♦  >  »  »  » 

component  of  >(T,Xq)  differed  from  the  corresponding  component  of  >(0;.,^^) 

by  more  than  10  The  Jacobi  constant  associated  with  the  solution  of  the 

equations  of  motion  remained  the  same  to  fifteen  significant  figu.es.  The 
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Jacobi  constants  associated  with  the  variational  equations  remained  the 

same  to  eleven  deciral  places.  The  equat ions (19)  were  satisfied  with 

-9 

residuals  at  most  10 

In  integrating  the  equations  of  motion  it  is  feasible  to  use  other 
methods  of  numerical  integration,  although  probably  more  machine  time 
would  be  involved.  The  power  series  method  becomes  relatively  more 
advantageous  when  Integrating  the  variational  equations.  The  solutions 
of  these  equations  are  highly  oscillatory  in  nature,  and  a  high  order 
numerical  integration  method  is  essential  if  an  extremely  small  step  size 
and  multiple  precision  in  adding  on  the  increments  is  not  to  be  required. 
Figure  1  shows  a  sketch  of  the  four  components  of  one  variational  equation 
solution.  The  points  showing  on  one  curve  are  thc»se  which  were  required 
to  be  calculated  by  the  power  series  method  using  16  terms  of  the  series. 

(In  actually  drawing  the  curves,  intermediate  points  were  also  calculated 
so  as  to  present  the  true  shape  for  illustration.)  Since  the  basic  orbital 
variables  occur  in  the  variational  equations,  the  step  size  required  for 
integration  of  the  variational  equations  is  then  the  step  size  for  the  whole 
problem  when  the  integrations  are  done  together  in  an  efficient  manner. 

In  many  types  of  problems,  loosely  approximate  solutions  of  variational 
equations  may  be  sufficient,  but  in  this  case  it  is  imperative  that  the  values 
of  the  four  variational  solutions  be  kno./r  ’"ight  at  the  end  of  the  period. 
Because  of  the  rapid  changes  in  the  four  solutio..s,  the  trace  of  the  latrix  of 
the  solutions  also  changes  rapidly,  and  the  caiculated  characteristic  roots 
would  not  be  correct  if  the  orbit  were  not  very  close  to  being  periodic  and  if 
the  variational  equations  were  not  solved  very  accurately.  In  order  to  be 
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absolutely  sure  that  the  accuracies  described  in  the  previous  paragraphs 
were  indeed  good  enough  to  permit  determination  of  the  characteristic  roots 
accurately,  we  chose  initial  conditions  which  were  approximations  to 
orbital  values  at  various  parts  of  the  orbit.  When  truly  periodic  orbits 
were  obtained  using  these  initial  conditions  as  first  guesses,  it  was  found 
that  the  calculated  characteristic  roots  were  the  same  to  eight  decimal 
places  as  t Key  were  when  a  start  was  made  at  a  different  point  on  the  orbit. 

In  accordance  with  the  criteria  given  in  Table  I,  the  results  in 
Table  VI  show  that  all  of  Rabe's  orbits  are  stable  except  the  last  one.  The 
stable  orbits  are  of  oval  shape.  The  one  unstable  orbit  is  of  horseshoe 
shape.  It  is  evident  from  looking  at  the  columns  listing  the  trace  and 
the  characteristic  roots  that  as  the  periods  of  the  various  orbits  increase, 
the  characteristic  roots  travel  around  on  the  unit  circle.  There  will  be 
points  of  "indifferent  stability"  as  indicated  in  Table  I  when  the  path  of 
the  roots  actually  hits  the  point  -1  or  -tl. 
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APPENDIX:  COMPUTER  PROGRAM  WRITE-UP 
Program  36  -  "Non-symmetr ic  orbits  with  variational  solutions" 

1.  General  Information 

A.  Purpose 

The  purpose  of  this  program  is  to  calculate  accurate 
periodic  orbits  in  the  plane  restricted  problem  of  thr# e 
bodies.  Four  independent  solutions  of  the  variational 
equations  may  also  be  computed  so  that  the  characteristic 
roots  associated  with  the  orbit  mav  be  obtained. 

B.  Restrictions 

1.  This  is  a  FORTRAN  IV  program  which  has  been  run  on 

an  IBM  7094  under  the  IBSYS  system  in  which  tape  5  is  the 
input  tape,  tape  6  is  the  output  tape  for  printing,  and 
tape  14  is  the  output  tape  for  punching. 

2.  All  input  variables  are  either  fixed  point  numbers  or 
double  precision  floating  point  numbers. 

3.  For  a  desired  orbit  the  exact  period  must  be  given. 
Approximate  initial  values  of  x,y,p^,  and  p^  are  also 
required;  arbitrary  starting  values  will  in  general  not  be 
good  enough. 

4.  If  one  asks  for  extreme  accuracy  in  finding  a  periodic 
orbit,  the  program  will  obtain  orbits  which  are  closer  and 
closer  to  being  periodic  for  a  while,  and  then  next  orbit 
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or  orbits  will  not  be  so  good  again.  The  reason  for 
this  is  given  on  page  34.  Thus  if  guesses  for  initial 
conditions  are  originally  extremely  good,  the  next 
supposedly  improved  guesses  will  not  be  as  good.  In 
this  regard,  the  values  of  the  FORTRAN  variable  C7  has 
the  effect  of  determining  whether  or  not  these  factors 
are  a  problem  in  a  particular  case. 

C.  Method 

1.  Mathematical  Method 

This  is  described  in  the  main  body  of  the  document. 
However,  as  stated  there,  some  of  the  actual  equations 
used  in  the  computer  program  are  slightly  different  from 
those  given  in  the  body  of  the  document.  There  are  4 
dependent  variables  in  equations  (3);  in  the  4  independent 
solutions  of  the  variational  equations  (9)  with  initial 
conditions  (18),  16  more  dependent  variables  occur.  In 
using  the  recurrent  power  series  method  for  solving  the 
20  equations,  we  here  introduce  13  auxiliary  variables  as 


follows  (instead  of  17  as  described  earlier): 
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(34) 


a  =  (x  + 

d  =  (x  +  u)y 

w  E  (x  +  u)^  + 

_ 1 _ 

®  "  (x  +  u)^  +  y^ 


(x  +  u  -  1)^  + 

A  5  1  -  _ 3(X  4-  U)^ 

(x  +  ;j)^  +  y^ 


B  =  1  -  — 3_Cx  -f  p  -  1) 

(x  +  p  -  1)^  +  y^ 

C  =  4-  p)v 

■  2  2 
(x  +  u)^  -♦.^y^ 

D  5  -f  u  -  l)v 

(X  +  p  -  1)^  +  y2 

R  Ml  -  U)[(x  +  p)2  + 

S  =  [(x  +  p  -  1)^  + 


L 


_ g-  r 

((K-^u)2  +  y2l^«L 


-  3(x-<-p)^ 
(x+p)^  +  y^- 


+ 


( (x+  p -1)^  + 


y2,3/2 


1 


_3(x4-p-1)^  1 

(x+p-1)^  +  y^J 


N 


3(l-p)(x-t-p)v  3p(x4-  p  -1)v 

[(x+p)^  +  ((x+p  -1)^  + 


If  these  definitions  are  written  in  the  form 
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(35) 


+  2  yX  +  vj 


d  ■  xy  +  u  y 
^  2 

w  ■  a  +  y 


g 

h 

A 

B 

C 

D 

R 

S 

L 

N 


1/w 

_ _ 

w  -  2x  +  (1  -  2u) 

1  -  3ga 

1  -  3ha  +  6hx  -  3(1  -  2ij)h 


gd 

hd  -  hy 
(1  -  u)g 
.  3/2 


3/2 


RA  +  SB 
3RC  +  3SD, 


it  is  seen  that  in  most  of  the  definitions  the  right-hand 
side  only  consists  of  linear  combinations  of  products  or 
quotients  of  at  most  two  of  the  other  variables.  In  the  case 
of  R  and  S,  the  equation  mav  be  differentiated  so  that 

2gR  *  3Rg 


(36) 


2hS  »  3Sh 
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which  will  also  be  satisfactory  for  the  recurrent  power 
series  method. 

In  all  we  have  33  variables  to  consider,  (20  desired 
dependent  variables  as  well  as  the  13  auxiliary  variables 
discussed  above).  Each  of  these  variables  is  considered 
as  a  power  series.  For  example. 


and 


■2 

n  »  1 
00 

■.? 


A  (At) 
n 


X  (At) 
n 


n-1 


n-1 


If  convenient,  the  coefficients  of  the  power  series  use 
the  same  symbol  as  the  dependent  variable  itself  but  are 
subscripted  with  an  n.  The  exceptions  to  this  show  up 
in  the  definitions 


P 


X 


p_,(4t) 


n-1 


2 

n»l 


,n(At) 


n-1 


6x^^^  -  5]  k^^^At)"'^  i  '  I, II, III, IV, 


,  (i)  V  (i),,,,n-l 

6y  ■  I  (At) 


i  -  I, II, III, IV, 
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(37.1) 

(37.2) 

(37.3) 

(37.4) 

(37.5) 


OD 


n 


i  -  I, II, II, IV, 


-  2  i  -  I, II, III, IV. 

^  n-1  " 

When  these  power  series  are  substituted  into 
equations  (35)  [or  (36)]  and  into  the  equations  of 
motion  (3)  and  in  the  variational  equations  (9),  and 
the  corresponding  powers  of  At  are  equated,  the 
following  recursion  formulas  are  obtained: 


nx  .  p  -f  y 
n+1  n  n 


n+1  ^n  n 


^  A  V2-j 


(37.6) 
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(37.7) 


(37.8) 


(37.9)  -  uy  .,  ^  X  V 


J  n+2-j 


(37.10)  .  -3^  «jV2-J 


(37.11)  =  -3(l-2i-)h  ,  -  3 


^  '"j  n+2-j 


(ha  ^  .  -  2h . X  . »  . ) 

j  n+2-j 


(37.12)  C  •  Za 


rn+2-j 


(37.13) 


Vl  ‘  -  N  Wj^ 


(37.14) 


Lp.,!  •  ^  V2-J^ 


(37.15) 


N  .,  «  3  A  (RjC^+2-J  ^jV2-J^ 


"Pnfl  “  '’n  ■  ^ 


n-H-j 


(37.16) 


(37.17)  “  -Pn  ^ 


(37.18) 


n+1  n  n 


i  -  I, II, III, IV 


(37.19) 


„,»> . .,»)  *  „(i). 

n+1  n  n 


1  -  I, II, III, IV 


(37.20) 


(37.21) 


„u»>  .  v<‘> 
n+1  n 


nv'iJ  -  -u<‘>  ^ 


1  -  I, II, III, IV 


+  ^  N 


1  -  I, II, III, IV. 


The  known  initial  conditions  are  the  quantities  y  v,,p 


i  -  I, II, III, IV.  The. 


der , 


(38.1) 


aj^  -  (Xj  +  u)‘ 


(38.2) 

(38.3) 


W,  -  a,  +  y, 


g,  -  1/Wj 


(38.4)  h 


1  ■  ^*1  +  1  -  2iJ 


(38.5) 


-  (1  -  u)r^ 
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(38.6) 

S  ■  u  h,  /hT 

1  1  1 

(38.7) 

^1  *  (Xj  +  u)yj 

(38.8) 

A,  -  1  -  3gja, 

(38.9) 

Bj  ■  1  -  3h^(aj  -  2xj  +  1  -  2ij) 

(38.  .:o) 

b  ‘  *1“] 

(38.11) 

b  *  •’I'b  -  ''!> 

(38.12) 

Lj^  ■  ^1^1  ^  ^1^1 

(38.13) 

Nj  -  3(R^C^  +  S^D^) 

Now  equations  (37)  may  be  used  in 

obtain  the  second  coefficients  In 

the  order  listed  to 

each  series,  etc.  When 

the  program  is  being  run  in  the  mode  in  which  cnly  the 
equations  of  motion  (and  not  the  variational  equations)  are 
being  solved,  equations  (37.9)  through  (37.15)  and 
equations  (37.18)  through  (37.21)  as  well  as  equations 
(38.7)  through  (38.13)  are  not  used  in  the  recurrer.ee 
process . 

For  the  mathematical  method  of  improving  the  "periodic" 
orbit  and  for  obtaining  the  characteristic  roots  associated 
with  the  periodic  orbit,  see  the  explanation  in  the  main  tody 


of  the  document . 
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2.  Coding  Method 

IBM  7094  FORTRAN  IV.  This  program  and  its  subroutines 
make  maximum  use  of  COMMON  storage.  We  have  desired  to  run 
many  similar  types  of  problems  without  changing  much  of  the 
program.  For  this  reason  there  are  many  COMMON  variables 
available  which  w^re  not  actually  used  in  this  version  of 
the  program.  These  show  up  as  being  undefined  in  the  listing 
of  what  each  COMMON  variable  represents. 

COMMON  Variables 

Note ;  A(1,I)  thru  A(9,I)  are  usually  as  defined  below. 

However,  at  the  end  of  the  orbit  under  the  control 
of  subroutine  ENDOR,  they  are  used  as  matrix  elements. 


A(1,I)  thru  A(12,I) 
A(13,I)  thru  A(24,I) 
A(25,I)  thru  A(33,I) 
A(34,I)  thru  A(49,I) 
A(50,I) 

B(l)  thru  B(12) 

B(13)  thru  B(20) 


^  i  ’  ^  1  ’  ^  1  ’  i  * 


J  I  I 


, II  „II  II  II 

k,  .tj  .Uj  ,Vj 


.III  „III  III  III  ,  IV  IV  IV  IV 

kj  .Ij  ,Vj  .  kj  .Uj  ,  aj.Wj.gj.h, 


^i  »  ^i  ’ 


®i’^i’^i’^i 


N, 


Ri^S^ 


^  I  I  A  I  I 
x.y.Px.P„*  .fip^.fip  , 

Ay  X  y 


^  II  ^  II  ^  II  .  II 
6x  ,(^y  ,6p^  ,cp^ 


^  III  .  Ill  ^  III  ^  III  ^  IV  ,  IV  ,  IV  ,  IV 
6x  ,6y  ,6p  ,6p  ,  ,6y  ,6p  ,fp 

Ay  X  y 


B(21)  thru  B(50) 


5b 


C(l)  thru  C(5)  x(0),y(0),p  (0),p  (0),  Jacobi  Constant  at  t«0. 

X  y 

C(6)  thru  C(50) 

E(l)  Trace  of  monodromy  matrix 

E(2)  Real  part  of  1st  characteristic  root 

E(3)  Imaginary  part  of  1st  characteristic  root 

E(4)  Real  part  of  2nd  characteristic  root 

E(5)  thru  E(50) 


F(I) 

GI(I)  i 

GII(I)  1/i 

S(I)  (At)^'^ 

x(I) 

y(I) 

z(I)  Erasable  storage.  (However,  at  the  end  of  an 

orbit,  z(I)  ^carried  to  subroutine  ORBIT 
through  subroutines  MATR,  GUES,  and  ENDOR) 

C1,C2 

C3  Parameter  used  to  help  determine  At.  (Making 

C3  smaller  improves  accuracy  attainable).  C3 
is  usually  chosen  between  .1  and  1.  and 
mostly  about  .25. 
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C4 

CS  Period  of  the  orbit 

C6 

C7  Accuracy  required  in  calling  orbit  "periodic". 

(This  is  the  allowable  difference  between  any  of  the 
4  initial  values  of  the  orbital  variables  and  the 
corresponding  value  at  the  end  of  a  period.) 

C8  thru  CIO  _ 

D  Max  At  which  subroutine  POWER  thinks  should  be 

used 

DT  At  which  is  actually  used 

DTMAX  Maximum  At  ever  to  be  used 

DSMAX 

El  thru  ElO 

G 


GMU 

y 

GMUl 

u 

-  1 

GMUC 

1 

-  2u 

H 
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I,II,J,K,K1  thru  KIO 

Erasable 

L 

M  -  1 

LI 

N  +  1 

L2 

N  +  2 

L3  thru  L8 

L9 

Usually  0.  If  1,  subroutine  CUES  prints  dump 
of  matrices  and  solutions. 

LIO 

Max.  number  of  times  subroutine  ORBIT  may  be 
called. 

M 

Number  of  terms  used  in  Taylor  series. 

Ml 

Case  number.  Later  this  is  used  by  subroutine 

ORBIT  to  tell  whether  it  is  through  (Ml  =  1) 

or  not  (Ml  *  0) .  j 

M2 

If  M2  *  4,  only  the  4  equations  of  motion  are 
solved. 

If  M2  «  20,  the  variational  equations  are 
solved  also. 

M3 

Number  of  integration  steps  per  printing  interval. 

M4  thru  M6 

M7 

(1,  there  is  to  be  no »  ^  ,  , 

If  )  f  intermediate  printing  of 

I2,  there  is  to  be  f  variational  quantities. 

M8 

M9 

Max.  number  of  times  subroutine  POWER  may  be  called 

MIO 
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N 

N1 

N2 


Running  variable  (N  *  1,M) 

If  N1  program  prints  | 

at  intermediate  points  on  o^'^'H. 

If  N2  variables  are  calculated. 

(2,  all  201 


N3 


A  counter  to  compare  with  M3 


N4 

N5 

N6 


If  N4 


{?}• 


subroutine  PRNT 


called  by 


subroutine  ORBIT. 


Ordinarily  N5  ■  0.  Subroutine  SETDT  sets  N5  *  1 
when  it  has  found  a  suitable  At  which  will  just 
end  an  orbit. 


N7 

N8 


If  N7 


there  is  now  no 
there  is  now 


intermediate  printing  of 
variational  quantities. 


N9  Max.  number  of  times  POWER  may  still  be  called. 

NIO 

p  Maxell 

Q.R 

T  independent  variable  (time) 

U.V.W 


xl  thru  xlO 
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yl  thru  ylO 

zl  thru  zlO  Erasable  storage 

II.  Running  Procedure 

A.  There  are  many  subroutines  used  by  this  program.  In  the 
order  that  they  are  called,  the  basic  four  are: 

RE36D  —  This  routine  reads  the  necessary  input  data  for  the 
problem  and  prints  it  back  out  again. 

SETUP  —  This  routine  initializes  the  necessary  quantities  to 
start  an  orbit.  It  prints  headings  and  calls  PRNT  to 
print  Initial  values  in  the  table. 

ORBIT  —  This  routine  runs  a  complete  orbit.  In  doing  so  it 
uses  the  subroutines: 

POWER  —  This  routine  calculates  the  series  coefficients 
and  computes  a  likely  size  for  D  =  At  . 

lllaX 

SETDT  —  This  routine  decides  just  what  At  to  actually 
use.  When  needed  it  calls  subroutine. 

TAYL  —  This  subroutine  calculates  the  dependent 
variables  for  a  given  At  by  the  Taylor 
series. 

SETPR  —  This  routine  decides  whether  any  printing  is 

to  be  done  at  the  particular  time  in  accordance 
with  input  parameter  information.  If  so,  it  tells 
the  program  ORBIT  to  call . 


IT 
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PRNT  —  This  routine  prints  the  orbital  parameters 
for  the  given  value  of  t. 

ENDOR  —  At  the  end  of  the  orbit,  this  routine  calculates 
and  prints  residuals  in  checking  equations  (19). 

It  c;-lculates  and  prints  the  trace  and  the 
characteristic  roots.  Then  it  calls  CUES  — 

This  routine  sets  up  equations  (12).  In  order  to 
solve  them  it  calls  subroutine  MATR.  Tht  suggested 
corrections  printed,  and 

the  new  values  of  C(I)  I  *  1,4  are  set  up. 

The  subroutine  ORBIT  then  decides  if  it  has  obtained  a  good 
enough  periodic  orbit.  It  transfers  this  information  to  the 
main  program  by  means  of  parameter  Ml. 

FANPR  —  This  subroutine  prints  all  necessary  information  (which 
might  be  saved  for  long  periods  of  time)  about  the  true 
periodic  orbit. 

The  main  program  also  punches  (actually  writes  output  tape  14) 

two  cards  which  give  Ij  ,  x  (0)  ,y  (0)  ,  Per  iod  ,  p  (0),p  (0)  in  the  format 

X  y 

required  as  input  to  the  program. 

B.  Data  input  cards. 

There  are  5  input  cards  required  for  each  case.  Any  number  of 
cases  may  be  run  one  right  after  the  other.  The  cards  are  as 
follows:  (When  a  quantity  is  given  as  "c. .  bi  rary"  it  means  that  this 
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input  quantity  will  never  actually  be  used  in  this  particular 
version  of  the  program.  Of  course,  if  one  wishes  to  always 
use  only  this  exact  version  of  the  program,  subroutine  RE36b 
could  be  easily  changed  so  as  to  read  only  the  quantities 
actually  needed.) 


Card 

1: 

Format 

(3D24.16) 

u 

x(0)  y(0) 

Card 

2: 

Format 

(3D24.16) 

Period 

Px^O)  Py(0) 

Card 

3: 

Format 

(1415) 

M  *  no . 

of  terms  in  Taylor  series 

M2  *  4  or  20  (20  means  variational 

equations  also  computed) 

M3  =  no.  of  integration  steps  per 
printing  interval 

M4  =  arbitrary 

MS  ■  arbitrary 

M6  “  arbitrary 

M7  ■  1  or  2  (2  means  there  is  to  be 

intermediate  printing  of  varia¬ 
tional  quantities) 

M8  “  arbitrary 

M9  *  max.  no.  of  times  subroutine  POWER 
may  be  called  (i.e.,  the  max.  no. 
of  integration  steps  you  wish  to 
allow  without  giving  up  on  ever 
getting  to  the  end  of  the  period). 

MIO  =  arbitrary 

LlO  =  Max.  no.  of  times  subroutine  ORBIT 
may  be  called  (i.e.,  how  manv 
complete  orbits  are  you  going  to 
allow;  if  LlO  =0  it  will  run 
one  orbit ) . 
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L9  *  0  or  1  (1  means  CUES  prints  a 
dump  of  matrices  and  solutions) 

L8  ■  arbitrary 

L7  ■  arbitrary 

Card  4:  Format  (7D10.?)  Cl  ■  arbitrary 

C2  «  arbitrary 

C3  *  parameter  used  to  help  determine 
At.  (Making  C3  smaller  tends 
to  cut  At).  C3  is  usually  chosen 
between  .1  and  1.  and  mostly 
about  .25. 

C4  *  arbitrary 

C6  ■  aibitrary 

C7  *  accuracy  required  in  called  orbit 
"periodic".  This  is  the  allowable 
difference  between  any  of  the  4 
initial  values  of  the  dependent 
variables  and  corresponding  values 
at  the  end  of  a  period. 

C8  *  arbitrary 

Card  5:  Format  (4D10.2) 

C9  “  arbitrary 
CIO  ■  arbitrary 

DTHAX  «  max.  value  of  At  you  would  ever 
want  it  to  use 

DSMAX  ■  arbitrary 

C.  Output 

Printed  output  depends  on  input  pa’-ameters  M3,  M7 ,  and  L9. 

Subroutine  RE36D  always  prints  the  input  data  it  has  read.  On 


I 
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a  new  page  will  then  come  a  table  giving  various  orbital 


prof iles . 

(If 

M3  is  very  large 

,  one  may 

only 

get  initial 

values  and 

values  at  the  end  of  a 

period) . 

The 

table  will 

generally 

be  in 

the  format 

T 

X 

y 

• 

X 

• 

y 

C 

6x' 

6y^ 

X 

6p^ 

V 

c‘ 

ix" 

6y” 

11 

op 

X 

.  II 

«Py 

c" 

6x'" 

A 

6y 

.  Ill 
“^Px 

<  HI 

6p 

y 

c’" 

,  IV 

OX 

.  IV 
oy 

A  IV 

*Px 

A  IV 
6p 

y 

■JV 

c 

T 

X 

y 

• 

X 

• 

y 

c 

6x' 

6y^ 

X 

6pl 

y 

c' 

«x" 

6y“ 

A  II 

6p 

X 

y 

c” 

«x”‘ 

A  m 

6y 

6p"' 

X 

^  III 

6p 

y 

c’” 

6x''' 

A  IV 

6y 

Sp'" 

X 

A  IV 

6p 

y 

c*'' 

etc . 

Of  course  if  M7  ■  1,  there  will  be  no  intermediate  printing 
of  the  variational  solutions. 

At  the  end  of  an  orbit  which  the  program  thinks  truly 
periodic  (or  even  if  it  has  been  unsuccessful  in  really  finding 
a  good  one),  the  six  variational  equation  checks  are  printed. 


i 

( 

/ 

i 


I 
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(These  are  the  residues  obtained  when  checking  in 
equations  (19),  and  they  should  all  be  very  close  to 
zero) . 

Next  are  printed  two  characteristic  roots  and  the  truce 
of  the  monodromy  matrix.  Then  are  printed  suggested 
changes  in  initial  values  as  obtained  as  a  solution  of 
equations  (32).  Then  the  new  initial  guesses  are  printed. 

If  the  program  has  been  successful  in  finding  an  accurate 
periodic  orbit,  one  will  usually  not  be  interested  in 
these  "suggested  changes  in  initial  values"  or  in  the 
"new  guesses". 

Fif.-’lly  the  program  prints  a  listing  of  important  orbital 
parameters  suitable  for  saving  as  a  page  in  a  book  of  orbits. 
If  L9  ■  I  ,  the  print  out  will  be  interspersed  with  matrix 
print  outs  by  subroutine  CUES.  These  are  not  labelled,  and 
the  easiest  way  to  see  what  they  mean  is  to  read  the  prog’'am 
listing. 

Punched  output  always  consists  of  two  cards  in  format 

(3D24.16)  giving  for  the  latest  orbit  u,  x(0),  y(0).  Period, 

p  ('^1.  p  (0). 

X  V 

D.  Results 

An  example  giving  the  print  out  for  a  particular  trojar.  orbit 
is  given  in  the  following  pages. 
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rRG36  OS/OT/RS 

EXTERNAL  FORNULX  NUNBER  -  SOURCE  STATEMENT  >  INTERNAL  FORMULA  NUNBERISI 


C  PROGRAM  36  NON-SVMMETRIC  ORBITS  MITH  VARIATIONS 

DOUBLE  precision  A,  B.  C.  Clt  C2.  C3.  CA.  CS.  C*.  CT,  C«t  CBt  ClOt 
10,  OT.  OTMAX.  OSMAX,  E,  El.  E2.  E3>  EA.  E5.  E6.  ET,  EB,  EB,  ElO. 
2F,  G,  Gl,  Gll,  GMU,  GMUl.  GMUC,  H,  P,  0.  R.  S.  T,  U.  V,  H,  XI.  X2. 
3X3,  XA.  X5.  X6.  XT,  X8.  X«,  XIO,  Tl,  T2,  V3.  VA.  VS.  V6,  VT,  TB, 

AY9,  YIO,  2.  21.  22.  23.  2A.  2S,  26,  2T.  2B.  29,  210 

DIMENSION  AI50.SOI,  BlSOi,  CISOI,  EI50I.  FI90I.  GIIMI.  GIIISOI. 
1X12000),  Y(2000l,  2(1001,  SISOI 
COMMON/SPR/X,  Y 

COMMON/OPR/A.  B,  C.  Cl.  C2.  C3.  CA,  CS,  CA.  CT.  CB,  C9,  CIO,  D, 
lOT,  OTMAX,  OSMAX.  E.  El,  E2.  E3,  EA,  E9.  E6.  ET.  EB,  E9.  ElO,  P. 

2G.  Gl,  Gll.  GMU,  GMUl,  GMUC,  H,  P,  Q,  R,  S,  T,  U.  V,  Ht  XI,  X2. 

3X3,  XA,  XS,  X6,  XT.  XB,  X9,  XIO.  Tl.  V2.  Y3,  YA,  TS.  YA,  YT,  YB, 

AY9,  YIO,  2,  21,  22.  23.  lA,  29,  2A,  2T.  2B.  29,  210 

COMMON/ 1  NT  S/1.  II.  J,  R,  Rl,  R2,  R3,  RA,  RS.  RA,  RT,  RB.  R9.  RIO. 
IL.  LI.  L2.  L3.  LA,  L5,  LA.  LT,  LB,  L9,  LIO.  M,  Ml,  N2.  M3,  NA.  NS, 
2N6,  MT.  MB.  M9,  MIO,  N,  Nl.  N2,  N3,  NA,  NS,  NA,  NT,  NB,  N9.  NIO 

1  FORMAT  (IHO  5023. IS) 

2  FORMAT  I12H0NEW  GUESSES  I 

11  FORMAT  (302A.16) 

12  format  IlHl) 

3  CALI  RE360 
k»1.0*30 

A  CALL  SETUP 
CALI  ORBIT 
IF  (Ml)  A. A, 20 

20  00  21  1«1.A 

21  2(1I-CII)«2II) 

WRITE  (A, 2) 

WRITE  16,1)  I2III.I-1.AI 

call  fanpr 


.B 

.10 

.12  .19  .lA 


WRITE  (A, 12) 

WRITF  I1A,11I  GNU,CI1I,CI2I,C5,CI9I,CIAI 

GO  TO  3 

END 


IB  .19 
20  .21  .22 
29 
2A 


BEGIN  ASSEMBLY  1A.020 


.IS 
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sumoutine  enoor 

DOUBLE  PRECISION  *•  B,  Ct  Cli  CZt  C)i  CR.  CSt  C6t  CT,  CB.  CBt  C10> 

10.  or,  OTMtR,  OSNRX.  E,  El,  E2.  E3,  E«,  E9,  E6,  ET,  EB,  EBf 

2E,  G,  Gl,  Gll,  GMU,  GNUl,  GMUC .  H,  P,  Q,  R,  S.  T,  U,  V,  M,  II,  12, 
313.  14,  15,  16,  IT,  18,  19,  110,  VI,  V2.  V3,  V4 ,  V5.  T6,  VT,  V8, 

4V9,  VIO,  {,  n.  12.  2  1,  24,  25,  26,  2T,  28,  29,  210 

DIMENSION  4150,501,  RI50I,  CI50I,  EI50I,  FI50I,  011501,  Glll50l, 
1II2000I,  VI2000I,  211001,  SI50I 
rONMON/SPR/l,  V 

COMMUN/UPR/A,  B,  C,  Cl,  C2.  C3,  C4,  C5,  C6,  CT,  CB.  C9,  CIO,  0, 
lOT,  OTPAI,  DSMAI,  E,  El,  E2,  E3,  E4,  E5,  E6,  ET,  EB,  E9,  ElO,  P, 

2G,  Gl,  Gll.  GMU,  GMUl,  GMUC.  H,  P,  0.  R,  S,  T,  U,  V,  H,  II,  12, 

313.  14,  15.  16.  IT.  18,  19.  110,  VI,  V2.  VI,  V4,  V5,  V6,  vf,  V8. 

4V9,  VIO,  2.  21.  22.  23,  24,  25,  26.  2T.  28,  29,  210 

COMMON/INTS/I .  II,  J.  I.  11,  12.  K3,  14,  R5,  K6.  RT,  R8,  IB,  RIO. 

11,  11,  12.  13.  14,  L5,  L6,  LT,  L8,  L9,  LIO,  M,  Ml,  M2,  M3,  M4,  N5, 
2M6.  MT,  MS,  M9,  MIO,  N,  Nl.  N2,  N3.  N4,  N5,  N6.  NT,  NB,  N9,  NIO 


1  FORM4TI33HOSII  VARIATIONAL  EQUATION  CHECRS  6012.51 

2  FORMAT  I IHOI 

3  F0RMATI112H0CHARACTERISTIC  EIPONENTS 

1  TRACE 

4  FORMATI 1H0202 1.14, 028. 14, 021. 14,028. 141 

NT  •  2  •> 

CALL  PRNT  ,2 

00  8  I >1,  3  ,1 

B2MI  >  BI5l*BII»13|,BI9l*BII»lTI-Bll3l>BII*5l*BllTt«BII*9t  ,4 

2121  •  2121-1.  ,8 

00  9  1-4,  5  ,T 

9  2111  •  BI6I*BI  1,1 1  l*BI  IOI*BI  l*15)-BI  14MBI  l«3l-BII8l*BI  l*Tt  ,B 

2151  ■  2151-1.  ,10 

2161  •  B(TI*R(16l*Bllll>BI20l-nil5l*BIBI-Bll91*BII2l  ,11 

WRITE  16.21  ,12 

WRITE  I6.1II2II),I<1,6I  ,14 

21  •  BI5|,BI10MBII51, 81201  ,19 

22  •  .25»ll-l.  ,20 

WRITE  16,31  ,21 

24  •  22,22*1.  ,23 

IF  1211  12,  30,  10  ,24 

10  IF  1221  20.  40.  12  ,25 

12  23  •  0SgHT(21*22l  ,26 

26  •  24-23  ,2T 

24  •  24*23  ,28 

GO  TO  45  ,29 

20  23  •  -21*22  ,30 

25  •  OSQRTI23I  ,31 

26  '■  24  ,,J2 

2T  .  -25  ,33 

GO  10  4T  ,34 

30  24  •  -1.  ,35 

26  •  -1.  ,36 

GO  TO  45  ,3T 

40  24  •  1.  ,38 

26  •  I.  ,39 


,9 


.9 


,13 

,19 


,22 


,16 


,IT 


45 

25  •  0. 

,40 

2T  •  0. 

,41 

4T 

Elll*21 

.42 

EI2I*24 

.43 

EI3I*25 

,44 

EI4I*26 

.45 

WRITE  16,41  24,  25,  26,  2T,  21 

.46 

CALL  GUES 

.49 

NT  •  MT 

.50 

RETURN 

.51 

END 

.52 

,18 


BEGIN  ASSEMBLV 


05.614 


67 


•.HTBOUTINf  r/INPH 

iiii'  f'if  purr  isiON  A,  p,  c,  ri,  c?,  cj,  ca,  ri.  C6.  ct,  cn,  ct,  cio. 

11),  IT,  UlMtx,  PSMAX.  E,  ri,  e;.  O.  f4.  15,  f6.  (7,  fS,  79,  710, 

if . I,  Ml,  i.Mu,  r.mii,  r.Muc,  h,  p,  o,  h,  5.  t.  u.  v.  Wt  7i,  x?, 

7X},  tA,  X5,  X6,  XT,  xa,  X9,  XlO,  Yl,  Y?,  YJ,  Y4,  Y5,  Y6,  Y7,  YX, 
4Y9,  VIO,  I,  ll,  12,  ii,  lA,  lb,  Ih,  n,  78.  79,  710 
OIPENSIUN  *150,501,  8(501,  C  ( 50 1 ,  7(501,  7  (501  ,  &M501,  (.11(501, 
1X(70C0),  YI^OOOI,  7(1001,  S(50l 


CUMPON/SPP/ X ,  Y 

COPPON/OPR/ »,  8,  C,  Cl,  C7,  C3,  C5,  C5,  C6,  CT,  C8.  C9,  CIO.  0, 
lOT,  OTP*X,  nSN*X,  7,  71,  77,  71,  75,  75.  F6.  77,  78.  79.  tlO,  F. 

7G»  r.l.  Oil.  OMU,  CMIIl,  GMUC,  H,  8,  0.  8.  S.  T,  U,  V,  M,  X|,  X7, 

IXl,  X*,  X5,  Xfc,  XT,  X8,  X9,  XlO,  Yl,  Y7,  YJ,  Y4 ,  Y5,  Y6,  YT,  Y8, 

5Y9.  Yic,  7,  71.  77.  71,  7*.  15,  76,  77,  78,  79,  710 

C0PM0N/INTS7I ,  II,  J,  K,  Kl,  K7,  R 1,  R6,  R5,  R6,  R7.  88,  89,  8|0, 
IL,  11.  l7.  11,  15,  15.  16,  17,  IB.  19,  LIO,  H,  PI,  N7.  Ml.  M4,  N5, 
7M6.  NT,  P8,  M9,  MIO,  N,  N| ,  N7,  Nl,  N5,  N5,  N6,  MT,  N8,  N9,  NlO 


1  F0RMAT(65H0 

SUN 

1  JllP(Tt«  I 

7  FORMAT (61M0 

CAR 

ITM  M0C6  1 

1  F0RMATt65M0 

EOU 

1*1  MASSES  1 

5  F0«M*T(51M0 

RM 

•  025 

1.161 

5  FORMAT ( 51M0 

JACOBI  CONSTANT 

■  025 

1.161 

6  FORM*  T ( 5 IMO 

R7RtOD 

>  025 

1.161 

7  F08MAT(56HOOM7  OF  TH7  CONJUGATF  COM817X  CH«8*C T 78 1  ST  1C  800TS  IS 
17075. 161 

8  70PMAT(56M0  TH7  TWO  8F«l  CHA88C T 78  I  ST  1C  BOOTS  *87 

17075. 161 

9  FURPATI56H0  TH7  t8*C7  N*S 


1075.161 

10  FORMAT  I 109H0  X  V 

I  8  SU8  X  8  SUB  V  I 

11  FORMATIlTHOINITIAl  VAIU7S  5075.161 

12  F08MATI17H0  FINAI.  VAI.UFS  5025.161 
I J  FORMAT ( IMOI 

15  F0RMATI16H0INITIAL  OX70T  ANO  OY/OT  VALU7S  M787  2025.161 

19  FORMATI75H1  81AN7  RCSTBICT70 

1  THRfF  ROOY  PBORIFMI 

MRITt (6,151 

IF(OMU-l. 0-01120, 20, 75 

20  W8ITF (6.1  I 
GO  TO  50 

25  IF(r.MU-. 017500178. 28.  10 
28  M8ITF (6,71 
GO  TO  50 

10  IF(GMU-. 500150, 15,50 
15  H8ITC (6,11 
50  W8ITF (6,5IGMU 


1  ,2 

5  ,5 

8  .9 

10 
11 

12  ,11 

15  ,19 


.16 


MRITt I6.5ICI5I 

,17 

.18 

.19 

MRITf (6,6IC5 

,20 

.21 

.22 

WRITE (6,111 

.21 

.25 

MR (Tt (6,101 

.25 

.26 

MRITt  (6.11  I  (Cdl.l'l 

.51 

.27 

.28 

.29 

WRITE  16,171 IBI 1  1,  1*1 

,51 

.12 

.11 

.15 

7l*C(7l*C(ll 

.17 

77-C(5l-C(l 1 

.18 

WRITE (6.15171,77 

.19 

.50 

.51 

IF  (M7-5I  55.70,55 

.52 

55  WRITE  16,111 

.51 

.55 

IF (E  (  1 1 150,55,50 

.59 

55  WRI TE  (6,8IE 17 1 ,E 1 5  1 

.56 

.57 

.58 

GO  TO  60 

.59 

50  MR| Tt 16, 7  IF  17 1 .E 1 11 

.50 

.51 

.52 

60  WRITE I6,9IE 1 1 1 

.91 

.95 

.55 

TO  RETURN 
7  NO 
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SUMOUTINC  cues 

DOUBLE  DECISION  ft.  B.  C.  Cl.  C2.  C3.  Cft.  CB.  C4.  C7.  CB.  C«.  Clu. 
10.  DT.  OTMftl.  DSMftX.  E.  El.  E2,  £3.  Eft.  E3.  Cft.  ET.  EB.  E«.  ElO. 

2E •  C.  Cl.  Oil.  CMU t  CftUl.  CMUC t  H.  ft.  0.  ft.  S.  T.  U.  h.  11.  12. 

313.  Ift.  IS.  Ift.  IT.  IB.  IS.  110.  Tl.  V2.  V3.  Tft.  TS.  Vft.  TT.  VB. 
ftTS.  TIO.  2.  21.  22.  23.  2ft.  23.  26.  27.  2B.  2S.  210 
dimension  ftlSO.301.  BI30).  riftOI.  E(SOI.  ftlSOI.  CIISOI.  CtItBO). 
1112000).  VI2000I.  21100).  SiSOI 
COMMON/ Sftft/1.  T 

COMMON/Dftft/ft.  B.  C.  Cl.  Cft.  C3.  Cft.  C3.  Cft.  CT.  CB.  CS.  CIO.  0. 
lOT.  OiMftx.  OSMxx,  E.  El.  E2.  E3.  Eft.  ES.  Eft.  E7.  EB.  ES.  ElO.  F. 
2C.  Cl.  Cl  I  .  CMu.  C**U1.  CMuC .  H.  ft.  Q.  ft.  S.  T.  U.  V.  W.  II.  12. 
313.  Ift.  IS.  Ift.  17,  IB.  IS.  110.  Tl.  T2.  T3.  Tft.  T3.  Tft.  TT.  YB. 
ftTS.  TIO.  2.  21.  22.  23.  2ft.  23.  2ft.  27,  19,  IS.  210 
COMMON/ 1  NTS/ 1 .  II.  J.  ft.  Rl.  12.  R3.  Rft.  R3.  Rft.  RT.  RB,  RS,  RIO, 
IL.  LI.  L2.  L3,  Lft.  L3.  Lft.  L7,  LB.  IS.  LIO,  M,  Ml,  M2,  M),  Mft,  MS, 
2Mft,  M7,  MB.  MS,  Mio,  N.  Nl,  N2,  N3.  Nft .  N3 ,  Nft,  NT,  N«,  NS,  NIO 

1  FOftMITI 3ftM0SUCCESTED  CHANCES  IN  INITIAL  VALUES  ft020.13l 

2  FOftMAT(33HOEOUATIONS  TO  BE  SOLVED  FOR  CHANCES  IN  INITIAL  VALUES  I 

3  F0ftMATIlH001S.12.ftD20.12l 
7  DO  a  l-l.  ft 

DO  B  J-2,  3 

R  ■ 

a  AII.JI  -  BIRI 
DO  S  l-l,  ft 

S  A(l,i*ll  ■  All.l*l)-1. 

DO  12  l-l.  ft 
12  Al l.fti  ■  Cl  I >-BI I  ) 

IF  (LSI  IB,  IB,  lft 
lft  MftITE  Ift. 21 

MftITE  lft.3IIIAI|  .J).Ja2.ftl.l«l.ftl 

IB  CALL  MATH 

IF  ILSI  23.  23,  20 

20  WRITE  lft.3l  I  IA|  I  ,  Jl , Ja2.ftl. I•1.SI 

23  WRITE  Ift, 1)12111. I>l.ftl 

21*0.  ^ 

00  30  lal.ft 
22>0A6SI2I I ) I 
IF  121-221  27.30.30 
27  2I-22 
30  CONTINUE 

IF  (21-ftl  30,30,35 
33  ft«.5«F 
21>F/21 
DO  ftO  l-l.ft 
ftO  2III*2III«21 

WRITE  Ift. II  1211). 1*1. ft) 

CO  lU  ftO 
30  F-21 
ftO  RETURN 
END 


.1 

.2 

.3 

.ft  .3  .ft 

,7 

.B  .S 

.10 

.11  .12 

.13 

.lft  .13 

•  lft  .17  ,1B  .IS  ,20  .21 

.22  .23 


.2ft 

,23 


.2ft 

,27 

.20 

.2S 

.30 

.31 

.32 

,3ft 

.33 

,35 

.3ft 

,37 

.30 

,3S 

,ftO 

.ftl 

.ft! 

•  ft3 

,ftft  ,ft3 

,ftft 

,ft7 

.ftB 

,ftS 

.30  .31 

.32  .33  . 

,57 
.30 
,5S 
.60 


3ft 


55 


3ft 
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SUBROUTINE  «*TB 

double  precision  a.  B.  C.  Cl.  C2.  C3.  C4.  CS.  C6.  CT,  CB.  C9,  CIO. 
ID.  DT.  DTMAX.  DSPAX.  E.  El.  E2.  E3.  E4.  E5.  E6.  E7.  EB.  E9.  ElO. 
2E.  G.  GI.  Gil.  GNU.  GMUl.  GHUC .  H.  P.  Q,  R,  S.  T.  U.  V.  k.  XI.  X2. 
3X3.  XA.  X5.  X6.  X7.  X8.  X9.  XIO.  VI.  V2.  V3,  V4.  V5,  V6.  V7.  V8, 

4V9.  YIO.  2.  21.  22.  23.  24.  2S.  26.  27,  28.  29,  210 

DIMENSION  AISO.IC).  B(SO),  CISC),  EISO).  FI50I,  G1I50I.  GI1I50). 
1X(2COO),  VI2000),  2(1001.  SISOI 
COMMON/SPR/X.  Y 

COMMON/DPR/A.  P.  r.  Cl,  C2,  C3,  C4,  C5,  C6,  C7,  C8.  C9,  CIO,  0. 
IDT.  OTPAX.  OSMAX.  E.  El,  E2.  E3.  E4,  ES.  E6.  E7.  E8,  E9.  ElO,  F. 
2G.  Gl.  Gll.  GMU.  GMUl.  GMUC .  H.  P.  Q.  R.  S.  T.  U.  V,  k,  XI.  X2, 
3X3.  X4.  XS.  X6.  X7.  X8.  X9.  XIO.  Yl.  V2.  V3.  V4,  VS.  V6.  V7.  V8, 

4Y9.  YIO.  2.  21.  22.  23.  24.  2S.  26.  27.  28.  29.  210 

COMMQN/INTS/I.  n.  J.  K.  Kl.  K2.  K3.  K4.  KS.  K6.  K7.  R8.  K9.  KIO. 
IL.  LI.  L2.  L3.  L4.  LS.  L6.  L7.  LB.  L9.  I  10.  M.  Ml,  M2,  M3,  M6,  MS, 
2Me.  MT.  ms.  M9.  MIO.  N.  Nl.  N2.  N3.  N4.  NS.  N6.  N7.  NB .  N9.  NIO 


RS-5 

.1 

6  KS  •  KS-1 

.2 

IF  IKS-1)  3S.  3S.  7 

.3 

7  21  •  0. 

.4 

K6»S-KS 

,5 

K2*K6'»1 

,  6 

KA  •  K2«l 

.7 

DO  IS  I«K6.4 

,8 

22  •  0ABS(A(I.K2)I 

.9 

IF  (22-211  IS.  IS.  13 

.10 

13  21  •  22 

.11 

K3  «  I 

.  12 

IS  CONTINUE 

.13 

.14 

DO  17  J»K*.  6 

.15 

17  A(K2*S.J)  •  A(K3. J l/A(K3.K2) 

.16 

.17 

DO  20  1>K6.4 

.18 

IF  (I-K3)  18.  20.  18 

.19 

18  DO  19  J>KA.  t 

.20 

A(I.J)  •  A( 1 . J)-A( I .k2I*A(K2*S. J) 

.21 

19  CONTINUE 

.22 

.23 

20  CONTINUE 

.24 

.25 

IF  IK6-K3)  24.  2B.  24 

,26 

24  DO  2S  J>ka,  6 

.27 

2S  AIK3.J)  >  4IK6.J) 

.28 

.29 

28  GO  TO  6 

.30 

3S  2141  •  AI4.6)/AI4.SI 

.31 

DO  40  I'l.  3 

.32 

K  •  10-1 

.33 

Kl  •  6-1 

,34 

2  IK-6)  •  A  IK. 6) 

.35 

DO  40  J>K1.  S 

.34 

40  2IK-6)  •  2IK-6I-AIK. J)»2( J-1) 

.37 

.38 

RETURN 

.40 

END 

.41 

,39 
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sumoutine  one  it 

DOUBLE  ODECISION  i.  B.  C.  Cl.  C2.  C3.  C4,  CS.  C6.  CT,  CS,  C«.  CIO. 
ID.  OT.  nSHAl,  E.  El.  E2.  E).  E«.  EB.  E*.  ET.  EB.  E9.  ElO. 

2f,  C,  Gl.  Gll.  OHU.  CMUl.  GMUC.  H.  B.  U.  B.  S.  T.  U.  V.  tt.  II.  I2. 
313.  *•,  15,  16,  1/,  18,  19,  HO,  TU  V2.  13,  »4,  n,  T6.  Y7,  V8, 

4V9.  vio.  it  11.  12,  13.  14.  13.  16.  IT,  18,  19,  110 

DIMENSION  6130,30).  B(30l,  CI30I.  EI30t.  EI30t.  Glt30l,  GIMSOl. 
1II2O00),  vdOOOl,  IIIOOI,  S(30l 
C0140N/SPR/I,  » 

COMi'ON/DP*/*.  8.  C.  Cl-  C2,  C>.  C4.  C3.  f4.  CT.  C8.  C9.  CIO.  0. 
IDT.  OTMII,  0SN61,  E,  El,  E2,  E3.  E4.  E3.  E6.  ET.  E8.  E9,  ElO,  f, 
2G,  Gl,  Gll,  G8U ,  G4U 1 ,  GNUC ,  H,  8,  0,  8,  S,  T,  U,  V,  W,  II,  12, 
313,  14,  13,  16,  IT,  IS,  19.  110,  VI,  T2.  V3,  V4.  V3,  T6.  V7,  V8, 

4V9,  VIO,  1.  11.  12.  13.  24.  13.  16.  IT.  IS.  19.  HO 

COMMON/ INTG/I,  II,  J,  K,  Kl,  K2,  K3.  R4,  13,  16.  IT,  IS,  R9,  KIO, 
IL.  LI.  L2,  L3.  L4,  l3.  L6.  IT,  LBi  L9,  LIO.  M.  *<1.  M2,  M3,  M4,  M3. 
2M6,  MT,  MS.  M4,  MIO,  N,  Ml  .  N2 .  N3.  N4.  M3,  N6,  NT.  NS,  N9,  NIO 


1  EORMiT  I IH02021. 14. 028. 14,021 . 14,028. 141 


3  call  'Owes 

.  1 

call  SETOT 

.2 

T-T*OT 

.3 

311  CALL  Sr.T8« 

,4 

IF  1 N4  1  7 .  7 ,  6 

.3 

6  CALL  PANT 

t  6 

7  IF  1 NSI  11,  11,  8 

,7 

8  IF  IN3I  3,  3.  12 

,8 

11  M2-4 

,9 

L10»0 

.  10 

12  IF  IM2-4I  13.14.13 

,11 

13  CALL  EMOOR 

,12 

14  IF  ILIOI  20,20,13 

.  13 

13  11>0. 

.  14 

DO  132  l>1.4 

.  13 

12-  OABSIC  n  1-3U  1  1 

.  16 

F  (11-121  130.132,132 

,1T 

130  11>12 

.18 

132  continue 

,19 

.20 

IF  Ill-CTI  20,20,133 

.21 

133  Ml«0 

.22 

DO  15 T 

.23 

137  CIII-KII  «CIII 

,24 

,23 

16  L10«L10-1 

,26 

aeiubn 

.2T 

20  Ml  .  1 

.28 

GO  TO  16 

,29 

.CUD 

.30 
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SURROUTINE  POWER 

DOUBLE  PRECISION  «.  B.  C.  Cl,  C2.  C>,  C4,  CS.  CB.  CT,  CB,  C«,  CIO, 
10,  OT,  OTMtx,  DSMIil,  E,  El,  E2,  E3,  E4.  ES,  E6,  ET,  EB,  E«,  ElO. 

«  G,  Gl,  Gil,  GNU,  GMUl,  CMUC ,  H,  P,  0,  R,  s,  T,  U,  V,  W,  Ml,  12, 
3X3,  X«,  IS.  X6.  XT.  IB.  19,  XIO,  Tl,  T2,  V3,  T9.  TS,  T«,  T7,  TB, 

9T9,  no.  2,  21,  22.  23.  29.  25.  26.  2T.  2B.  29.  210 

DIMENSION  Also, SOI,  BISOI,  C I  SC  I ,  EISOI.  EISOI.  GII50I.  GlltSOI, 
11(20001.  YI2000).  21100),  SISOI 
COMMON/SPR/>.  V 

c  nM«()N/DPR/».  «.  C,  Cl.  C2.  C3.  C9.  CS,  C6.  CT.  CB*  t9,  CIO,  D, 
lOT,  DTmk,  OSMII.  E.  El.  E2.  E3.  C9,  ES,  Eb.  ET.  EB.  E9.  ElO.  f, 

2(j  t  Gl.  Gll.  GHU .  GMUl .  GMuC .  H.  P.  0.  R.  S.  T.  U.  V.  W.  XI.  12. 
3X3.  19,  IS,  16,  XT,  IB.  19,  XIO,  Yl,  T2.  V3.  Y9,  TS,  T6,  TT,  TB, 

9Y9,  vio.  2.  21,  22.  23.  29,  2S.  2b.  2T.  2B.  29.  210 

COMHON/INfS/l  ,  II,  J,  K,  Kl,  R2,  K3,  R9.  XS,  Rb,  RT.  RB.  R9,  RIO, 
IL.  LI,  L2.  L3.  19,  LS,  Lb.  LT.  LB.  19.  LIO.  H.  Ml.  N2.  M3.  M9.  MS. 
2M6.  MT,  MB,  H9,  Mio,  N,  Nl,  N2.  N3.  N9,  NS,  Nb.  NT,  NB.  N9,  NIO 


DO  10  I*).  M2  ,I 

10  X(  1 .1  I  •  Bl 1 1  ,2.3 

2  1  «  BIIUGMU  .9 

bl?l.  1  I  •  21*21  .5 

xi22.n  >  XI2I  .n*BI2l«BI2l  .6 

bl23,  1  I  •  l./XI22,l I  . I 

*129.11  •  l./l*(22. 11*1. -21-211  .B 

*I7S.!)  •  -GHU1**I23. 1  l•OSO*T 1*1  23,1  n  ,9 

*(2b.l)  •  GMU**I29, 1 l•0SCRTI*(29,  1  I  I  ,10 

*IS0.1I  •  *I2S.1 1**126,1 1  .11 

GO  TO  120.19),  N2  .12 

!<>  *  (2T.  1  I  •  21*B(2I  ,13 

*128,11  >  1.-3.**I23.1)**I21.1I  .19 

*129,11  •  1. -3. •*(  29,  1  )• (1(21 . 1  I  *1.-21-211  , IS 

*1 30,  1  I  •  *123, 1 I•*I2T, 1  I  ,16 

*131.11  •  *(29.1I*I*I2T,1I-B(23l  .IT 

*(32.11  •  *(2S-<<**(2B.1I**(2*.1I*|I29,1I  ,1B 

*(33,11  • (*(2S,1 )•*( 30, 1  I **(26, 1  I**!  31,1 ) 1*3.00  ,19 

20  00  SO  N>1,  L  ,20 

LI  -  N*1  ,21 

L2  •  N«’  ,22 

21  ■  0.  ,23 

22  •  0.  .29 

23  •  0.  .2S 

*(1.N*1I  •  GII(NI*I*I3,N)**(2.NI)  ,26 

*(2.N*1I  •  GII(NI*(*I9,NI-*(1.N) I  ,27 

00  29  J*t,  LI  ,2B 

*  «  L2-J  ,29 

21  ■  21**(1 .J)**( l.R)  ,30 

22  >  22**(1,JI**(2.RI  ,31 

29  23  -  23**(2.J1**(2.RI  ,32  .33 

*(?1,L1I  *  2.*Gmu**I 1 ,L1  )*21  ,39 

*122. Lll  ■  *I21,L1)*23  ,  3S 

29  ■  0.  .36 

25  •  0.  ,3T 

26  ■  0.  .3B 
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SUBROUTINE  POWER 

DOUBLE  PRECISION  R,  B,  C.  Cl.  C2.  C3.  C4.  CS.  C6.  CT,  CS,  C«.  CIO. 
ID.  OT.  OTNRX,  OSNRI.  E.  El.  E2.  E3,  ER.  ES.  E6.  ET.  Et.  E9.  ElO. 
2F.  C.  Cl.  Gll.  GHU.  GMUl,  GMUC.  H.  P.  0.  R.  S.  T.  U.  V.  W.  XI.  X2. 
3X3.  XR.  XS.  Xb.  XT.  X8.  X9.  XIO.  VI.  V2.  V3.  Tb.  VS.  V6.  VT.  VB. 

bV9.  VIO.  2.  21.  22.  23.  2b,  25.  2b,  2T.  2B,  29.  210 

DIMENSION  R(SO.SO),  BISOI,  CI50I,  EtSOi,  PISOI,  CIISOI,  GlltSOI, 
1X(?000).  VI2000),  21100).  SISOI 
COMMON/SPR/X.  V 

CnMMON/OPR/R.  M.  C.  Cl.  C2.  C3.  Cb.  CS.  Cb.  CT.  C6.  CIO.  D. 
lOT.  DTPbX.  OSHAX.  E.  El.  E2.  E3.  Cb.  ES.  Eb.  ET,  E8,  E9.  ElO.  F. 
2G.  G1 .  Gll.  GNU.  GMUl.  GMUC.  H.  P.  0.  R.  S.  T.  U.  V,  M.  XI.  X2. 
3X3.  Xb.  XS.  Xb.  XT.  XS.  X9.  XIO.  VI.  V2.  V3.  Vb.  VS.  Vb.  VT,  VS, 

bV9,  VIO.  2.  21.  22.  23.  2b.  2S.  2b.  2T.  28.  29.  210 

COMMON/INTS/I.  II.  J.  K.  Kl.  R2.  K3,  Rb.  KS.  Kb.  KT.  KB.  K9.  RIO. 
IL.  LI.  L2.  L3.  Lb.  LS.  Lb.  LT.  L8.  L9.  LIO.  M.  Ml.  H2,  N3.  Mb.  MS. 
2Mb.  M7.  MB.  M9.  MIO.  N.  Nl.  N2.  N3.  Nb.  NS.  Nb.  NT,  NB,  N9.  NIO 


DO  10  1-1.  M2  ,1 

10  X( 1.1 )  •  B( 1  I  .2.3 

21  •  Blll-GMU  .b 

*121.11  -  21*21  .5 

*127.11  •  *(21.n*B(2l*B(2>  .b 

*123. 1  I  •  l./*(22.1)  .  I 

*l2b.l)  -  1.21*122. 11*1. -21-211  .B 

*(7S.1>  •  -GMU1**(23.1)*DS0RTI*I23.1II  .9 

*126.1)  •  GMU**l2b.l)*0SCRTI*l2b,l))  ,10 

*IS0,1)  •  *l2S,l)**l2b,l)  .11 

GO  TO  (20.1bl .  N2  .12 

lb  *I2T.1)  «  21*BI2)  .13 

*128.1)  •  1.-3.**I23.1)**I21,1)  .lb 

*129.1)  •  l.-3.**12b. 1)*1*I21.1)*1.-21-21)  ,15 

*130.11  «  *I23.1)**12T.1 )  .16 

*131.1)  •  *12b.ll*l*12T.l)-BI2U  .IT 

*132.1)  •  *125.1)**128.1)**12b.l)*|l29.1l  .IS 

*133.1)  •1*I2S.1)**130,1 )**12b.l )**I31.1I )*3.D0  ,19 

20  DO  SO  N-1.  I  .20 

LI  •  N*1  ,21 

L2  •  N*2  .22 

21  ■  0.  .23 

22  •  0.  ,2b 

23  •  0.  .25 

*11,N*1)  •  GI I  IN)*I*I3,N)**12.N)  )  ,2* 

*12.N*1)  •  GlllN)*l*lb,N)-*ll.N)l  ,2T 

DO  2b  J-l,  LI  ,28 

R  •  L2-J  ,29 

21  •  21**11. Jl**l l.KI  ,30 

22  •  22**11. J)**12,RI  .31 

2b  23  •  23**12. J)**12.R)  ,32  ,33 

*121, LI)  •  2.*GMU**I1,L1 )*21  ,3b 

*122, LI)  •  *I21,L1I*23  ,SS 

2b  ■  0.  .3* 

25  •  0.  .3T 

2b  •  0.  ,38 
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SUBKOUTINE  raNT 

DOUBLE  DECISION  A,  B,  C.  Cl.  C2.  C).  C«.  C5.  C6.  CT.  CB.  CB.  CIO. 
10.  DT.  OIMAX.  OSMAX,  E.  El.  E2.  El.  EA.  EB.  E*.  ET.  EB.  E«.  ElO. 

2E .  G.  Cl.  Gll.  G^U.  GWUl.  GMUC.  H.  E.  C.  B.  S.  T.  U.  V.  W.  XI .  X2f 
1(3.  X4,  xs.  x«.  rr,  xs.  x«.  xio.  ti,  t2.  vi.  V4.  tb.  t*.  v7,  tb. 

«V4.  ''10.  2.  21.  22.  21.  24.  2B.  26.  27.  2B.  2«.  210 

OIHENSIUN  6IB0.B0I.  BIBOI.  CIBOI.  EIBOI.  ElBOt.  CMBOl.  6IIIB0I. 
1X120001.  T(2000).  2(1001.  SIBOl 
CO»»»»ON/SPX/x,  V 

COXMON/DPK/6,  B.  C.  Cl.  C2.  Cl.  C6.  CB.  C6.  C7.  CB.  CB.  CIO.  D. 
lOT.  OTPXX.  OSMXX,  E  El.  E2.  El.  E4.  ES.  E6.  E7.  EB.  EB.  ElO.  E. 
2G.  Gl.  GM.  GMU.  G(xul.  GBUC .  H.  E.  Q.  X.  S.  7.  U.  V.  W.  XI.  X2. 
1X3,  X4,  X3,  X6.  XT.  XB;  XB.  XIO.  71.  72.  71.  74.  7B.  76.  77.  7B. 

479,  710.  2.  21.  22.  23.  24.  2B.  26.  27.  2B.  2B.  210 

C0«N0N/INTS2I  .  II.  J.  K.  Kl.  K2,  Kl.  K4.  KB.  K6.  K7.  KB.  KB.  KIO. 
U.  LI.  L2.  L3.  L4,  lb.  L6.  LT.  LB.  L9.  LIO.  M.  Ml.  N2.  Ml.  N4.  MB. 
2M6,  mT,  mb,  M9.  MlO.  H,  Nl.  N2.  Ml.  M4.  MB.  N6.  N7.  MB.  MB.  NIO 

1  EOBMXTI1H001B.B.4D20.  12. 02 3. IB  I 

2  FOXHATdH  4022. 14.021. IBl 
N4  •  0 

22  •  6(21 ♦B' 31 

23  •  6(41-611 
21  •  6(2)*6(2i 

24  •  (6(1 l*GMUI«<2*21 

25  >  (6.11  GMUl l••2«21 

26  >  DSJRT(24) 

27  «  OBOI)T(25l 

29  .  GPU*(2./iT*2B»-GMUl«(2./26«24l-22«22-21*21 
MMITE  (6.1)T,  B(11.  B(2I.  22.  21.  29 
Ll  •  I3*l 
X(L3I  •  6(11 
7(131  ■  6(21 
GO  TO  120. M.  N7 
6  24  •  l./(24*26l 

25  •  l./(25«27l 

21  •  CEU4/4-GMU1»24 

26  •  6(4l-21«6(li*CMU«GMUl«(24-25l 

27  •  -6(ll-2l*6(2t 
00  a  l«l.  11.  4 

2a  •  -TB^BI  l*4l-27«B( I«BI*22«B(|46I*21«6II«7I 
a  NKITE  (6.216(1*41.  B(I*5I.  B(I«6I.  611*71.  It 
20  RETUKM 
END 


.7 

.6 

.B 

.10  .11  .12 

.11 

.14 

.IB 
.16 
.17 
•  IB 
.19 
.20 
.21 
.22 
.21 

.24  .29  ,26 

.26 

.29 


.27 


BEGIN  6SSEMB17  OB. 619 
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SUOROUTINF  PC  360 

ooilfllf  PHrCISION  *,  B,  C,  Cl,  C?,  C3,  C6,  C5,  C6,  CT,  C8,  €■»,  CIO, 
10,  OT,  OtBAX,  OSH**,  E,  El,  E2,  CS,  t6,  €5,  C6,  ET,  FB,  C9,  EIC, 

2F,  r.,  Ri,  r.ii,  r.Bu,  GMul,  r.Muc,  h,  p,  o,  b,  $,  T,  ii,  v,  w,  xi,  x2, 

3K3,  »*,  «?,  «6,  X7,  >8,  *■»,  *10,  VI,  V2,  Y3,  VB,  V3,  V6,  V7,  V8, 

♦  V<J,  VIO,  /,  /I,  72,  73,  74,  75,  76,  77,  78.  79,  710 
OIMINSIUN  »(50,S0),  B(50t,  C  I  50  ( ,  E  I  50  » ,  F(50».  01150).  0,11(50). 
1*(7000),  Y(7000).  7(100).  5(50) 

COHMCin/spb/X,  Y 

CtlMMON/OP*/*,  8,  C,  Cl,  C2.  C3,  C4,  C5,  C6,  CT,  CB,  C9,  CIO,  0, 


lOT,  OTPAX,  OSMAX,  E.  El,  72,  C3,  F4, 

E5, 

E6.  ET, 

ER. 

€9, 

E  10 

,  F, 

2G,  Gl  ,  Gll,  GM(|,  GMIIl,  GMUC ,  H,  P,  0 

,  R , 

S,  I,  U 

,  V, 

Wo 

XI. 

X2. 

3X3,  X4,  15,  X6,  XT,  XB,  X9,  XIO,  Yl, 

Y2, 

Y3,  Y4, 

Y5, 

V6, 

YT, 

YB, 

4Y9,  YIO,  7,  7’,  72.  73.  74,  75,  76, 

IT, 

78,  79, 

710 

COMMON/INTS/I .  II,  J,  R,  RI,  R2.  R3. 

X4. 

R5,  R6, 

XT, 

K0O 

X9, 

RIO. 

!(.,  11,  12,  L3,  14,  L5,  16,  LT,  LB,  19,  LIO,  M,  Ml,  M2.  M3,  M4 ,  M5, 
2M6,  M7,  MR,  M9,  MIO,  N,  Nl,  N2.  N3,  N4,  N5,  N6,  NT,  NB,  N9,  NIO 

1  F0PM»T(14I5) 

2  FO«IMAT17ni0.2) 

3  FOPMAT ( 7014,6 ) 

4  FORMAT  (2023.16) 

5  format  (3024.16) 


READ  (5,5)  GMll,C(ll,C(2l,C5,C(3l,C(4) 

.1 

.2 

,3 

REAO  1 5, 1 IM,M2,M3,M4,  M5,  M6 ,  MT,  MB,  M9,  M|0.  LIO,  L9,  L8,  LT 

.4 

.5 

,  6 

,T 

.0 

.9 

.10 

.11 

.12 

.13 

.14 

.15 

.16 

.IT 

.18 

.19 

READ  I5,2IC1.  C2,  C3,  C4,  C6,  CT,  CB,  C9,  CIO,  OTMAX,  OSMAX 

.20 

.21 

.22 

WRITE  (6,5)  GMII.CI  1  )  ,CI2I  ,C5,CI  3I,CI  4) 

.23 

.24 

.25 

WRITE  (6, DM, M2. M3, M4,  M5,  M6.  MT,  MB.  M9,  MIO,  LIO,  L9,  LB,  LT 

.26 

.2T 

.28 

WRITE  (6,3IC1,  C2,  C3,  C4,  C6,  CT,  CB,  C9,  CIO,  OTMAX,  OSMAX 

.29 

.30 

.31 

return 

.32 

ENO 

.33 

BEGIN  assembly  12.B79 


SUBROOTINE  SETOT 

OO(JB)  E  PRECISION  A,  B,  C,  Cl,  C2,  C3,  C4,  C5,  C6,  C7,  C8,  C9,  CIO, 

10.  OT,  OTMAX,  OSMAX,  E.  El,  E2,  E3,  E4,  E5,  E6.  ET,  EB,  E9,  ElO, 

2F ,  G,  Gl,  Gll,  GMU,  GMUl,  GNUC ,  H,  P,  0,  R,  S,  T,  U,  W,  M,  XI,  X2, 
3X3,  X4,  X5,  *6,  XT,  X8,  X9,  XIO,  Yl,  Y2,  Y3,  Y4,  Y5,  V6,  YT,  Y8, 

4Y9,  VIO,  7,  71,  72,  73,  74,  75,  76,  IT,  78,  79,  710 

01  ME  NS  I  ON  Also, 50),  8(50),  C(50),  E(50l,  F ( 50  I ,  Gll 501,  Gll 1 50  I, 
1*12000),  Y(2000),  71100),  S(50)  oiii:,ui, 

COMMON7SPR7X,  Y 

COMMON/OPR/A,  B,  C,  CU  C2,  C3,  C4,  C5,  C6,  CT,  C8,  C9,  CIO,  0, 
IDT,  OTMAX,  OSMAX,  E,  El,  F2.  E3,  E4,  E5,  E6,  ET,  E8.  C9,  ElO,  F, 
2G,  Gl,  Gll,  GMII,  GMUl,  r.MUC,  M,  P,  Q,  R,  s,  T,  U,  V,  M,  XI,  X2, 

3X3,  X4,  *5,  X6,  XT,  XB.  X9,  XIO,  Yl,  Y2,  Y3,  Y4,  Y5,  Y6,  YT,  Y8, 

4Y9,  VIO,  7,  71,  72,  73,  74,  75,  76,  IT,  78,  79,  710 

COMMON/ INTS/ I,  II,  J,  K,  RI,  R2,  K3,  K4,  K5,  R6,  XT,  K8,  R9,  RIO, 

11.  LI,  L2,  L3,  L4,  15,  16,  IT,  LB,  19,  110,  M,  Ml,  M2,  M3,  M4,  M5, 

rSt  7’’  NT,  NN,  N9,  NIO 

IF  (TF0“C5I  2,  6,  6 

2  OT  •  0 

3  CALL  TAVL 
RETURN 

6  OT  •  C5-T 
N5  •  1 
GO  TO  3 
END 


BEGIN  assembly 


03.862 
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suAROtirixr  sum 

ooum  PBirniON  •.  e,  Cf  ci>  12,  o,  c«.  c^i  ce,  cr.  rt.  C9,  ciOt 
m,  PI,  niMAi,  osHBi,  r.  n.  f2,  ft.  »«•  («•  22,  f,  (9.  no, 

22,  01,  r.li,  r.Mu,  r.Mul,  owk:  •  h,  p,  0,  a,  S,  r>  u,  v,  ■,  ii, 

ID,  I«,  I^,  16,  IF,  la,  19,  no,  TI,  yl,  VI,  V9,  t9,  Vk,  vr,  yp, 
*v9,  Tio,  /,  ii,  II.  IS,  f*,  IS.  IP,  It,  It.  19,  no 

oinNsiim  Al^o,90l,  aikoi,  cnoi.  moi,  vnoi,  oM^ol,  ctmoi, 
imoooi,  vnooci  /iiooi,  snoi 

co"«ON^SP»/i,  y 

coNMnN/naa/t,  a,  t,  ti,  t2,  tv,  c*,  IS,  Ck,  CT,  ca,  to,  tio,  n. 


lot. 

piMki ,  oiMki ,  r ,  F 1 

,  2  2,  11.  Fk, 

f  k. 

ca.  CT, 

ca. 

C9, 

1  10 

.  2, 

2G. 

Gl  ,  Gll,  GMII,  GMUl  , 

f.MOt,  M,  a,  0 

.  a. 

S,  T,  0 

.  V. 

9, 

■1. 

12. 

Ill, 

■k,  Ik,  16.  IF,  IR, 

19,  no,  TI, 

T2, 

»J,  Tk, 

Tk. 

Tk, 

TF, 

TI, 

kT9, 

TIO,  2 .  21.  12,  IS, 

2k,  IS.  IP, 

IT. 

20.  29, 

no 

C09M09/ 1911/ 1 ,  11,  J.  a 

,  ai,  a2.  ai. 

ak. 

ak,  Ik. 

a?. 

aa. 

a9. 

BIO 

U,  II,  l/,  O,  (4.,  Ik,  Ik,  IF,  la,  19,  no,  M,  Ml,  M2,  M),  H«,  Hk, 
2Mk,  MF,  MM,  «9,  aio.  N,  Ml,  92,  9).  9k,  9k.  9k,  9T,  90,  99,  9|0 
9)  •  91*1 

IF  I9I-9II  2,  k,  k 
2  9k  •  0 
1  It  l<W9 
k  9k  •  0 

9k  •  I 

GO  10  I 
f  90 

•EGIN  kSSEMRlV  Ok.Rhk 


^uaaouiiNf  SfTua 

ooiiBi.(  aRfciki09  «,  a,  c,  ci,  C2,  ci,  Ck,  ck,  Ck,  cr,  ca,  C9,  cio, 
10.  or.  OFMki.  osMki.  2.  n-  rr.  rv,  ik,  (k,  Ck,  cr,  la,  19,  Cio, 
22,  G,  Gl,  GM,  GMU,  GMU|,  GMUCt  H,  P,  Q,  a,  Sj  I,  U,  IF,  9,  II,  12, 
in,  Ik,  Ik.  16,  II.  la.  19,  no,  ti,  t2.  ti.  Tk.  Tk,  to,  ti.  ta, 
9T9,  tio,  I,  21,  12,  ti,  IP,  IS,  2k,  21.  2a,  29,  2tO 
OIMC9SI09  aiko.koi,  aikoi.  cikoi,  cikoi,  rikoi.  ciikoi,  ciiikoi. 

III7000I,  T(2000l.  211001,  SIkOI 
C0MMn9/kaaFi,  t 

coMM09/opa/a,  a,  c,  ci,  C2,  ci,  Ck,  ck.  ca,  cT,  ca,  C9.  cio,  o, 
lor,  oraai,  okMki,  c,  ci.  C2,  fi.  Ck.  ck,  ca.  cr,  ca.  C9.  no,  2, 

20,  Gl,  Gll,  GHU,  GMOI ,  GMUC,  H,  P,  0,  a,  S,  1.  U.  V.  M,  l|,  12, 

ii>.  Ik,  Ik,  la.  IT,  la,  19,  no,  ti,  v2,  ti,  Tk,  Tk,  va.  ti,  ti, 

kT9.  TIO,  2,  21,  22.  21.  2k.  2k,  26,  21,  20,  29,  210 

009909/1911/1.  II,  J,  a.  ai.  12,  a),  ak,  ak.  aa.  ar.  aa.  I9.  no, 
a,  II,  12.  LI,  Ik,  Ik,  Ik,  ir.  ca,  i9,  no,  h,  mi.  m2,  mi,  Mk,  Mk. 
296.  Ml,  Ma.  99,  Mio,  9,  Ml,  92,  Ml,  Nk,  9k,  96,  MT,  90,  99,  910 
2  COaMATI 1 1 2H|  I  I  T 

1  OI/OT  OT/Of  C  I 


T  •  0. 

.1 

9T  .  9? 

.2 

1  •  9-1 

.1 

L)  •  0 

,k 

9)  •  0 

,1 

Mk  •  0 

,6 

99  «  99 

.T 

Sill  •  1. 

.a 

GMtII  •  GMU-I. 

,9 

GMUC  •  l.-GMU-CMU 

.10 

DO  T  !•!,  k 

.11 

T 

a  1 1 1  •  cm 

.12 

,11 

00  8  m.  M 

.  Ik 

GUM  ■  1 

.11 

a 

GIMII  •  l./GMIl 

.16 

.IT 

92  •  1 

.la 

IT  192-kl  20,  20,  9 

.19 

9 

92  •  2 

.20 

Ik 

no  Ik  1*6,  19 

,21 

Ik 

81  1  1  •  0. 

.22 

,21 

DU  Ik  1*1,  20,  k 

.26 

16 

8111  •  1. 

.29 

•  26 

20 

MR  IFF  16,21 

.2? 

.21 

cm  RRMT 

,29 

Clkl-29 

.10 

RCTURM 

.11 

(Mn 

.12 

acciM 

kSSCMBlT  Ok. 621 
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SUB«OUTINf  T»YL 

DOUBLE  PRECISION  «,  6.  C.  Cl.  C2,  C3.  C4.  CS.  C6.  Cl,  CP.  C9,  CIO. 

ID.  OT,  OTMRX.  DSHRX,  E.  El.  E2.  E3.  EN.  EB,  E6,  ET,  E8.  E9,  ElO. 

2P.  C,  Cl.  on.  CHU.  GNUl,  CMUC,  H,  P,  0,  R.  S.  T,  U,  V,  R,  XI.  X2, 

3X3,  X9,  X5,  X6.  X7,  XP,  X9,  XIC,  Yl,  Y2,  Y3.  YL,  YS,  Y6.  Y7,  Y8, 

9Y9,  Yio.  I,  21.  22.  23,  29.  25.  26,  27,  28,  29.  210 

DIMENSION  A(5C',50).  B(50l.  C(50).  E  I  50 ) .  F  (  SO )  ,  GIISO).  GIKSC). 

1X(2000).  Y(2C0C),  211001,  SISC) 

COMmON/SPR/X,  V 

COMMON/DPR/*,  B,  C,  Cl.  C2,  C3.  C*.  C5,  C6,  C7,  C8,  C9,  CIO,  D, 
lOT,  DTM4X,  DSMXX,  E.  El.  E2,  E3.  E*.  E5.  66,  E7,  E8,  E9,  ElO.  F, 

2G.  GI.  Gll.  GP-U.  GMUl,  GMUC .  M.  P.  Q.  R.  S.  7.  U.  V.  R.  XI.  X2. 

3X3,  X9,  X5.  X6,  X7,  X8,  X9,  XIO,  Yl,  Y2,  Y3,  Y9,  Y5,  Y6.  Y7,  YS, 

*Y9,  YIO,  2,  21,  22,  23,  29,  25,  26.  27,  28,  29,  210 

COMMON/ ISIS/ I ,  II,  J,  X,  XI,  X2,  X3.  X9,  X5,  X6.  X7,  KB,  X9,  XIO, 

IL,  LI,  L2,  L3,  L9,  L5,  L6,  l7,  L8,  L9,  LIO,  M,  Ml,  M2,  «3,  M9,  m5, 

2M6,  M7,  mb,  MR,  MIO,  N,  Nl.  N2.  N3,  N9,  N5,  N6 ,  N7,  N8,  N9 ,  NIO 

,3 


.8  .9 


BEGIN  ASSEMBLY  03.869 


DU  2  1*2,  M  .1 

2  SI n  •  SI l-l)«DT  .2 

DO  8  I«l,  MZ  ,9 

Bl I  I  •  All , 1 1  .5 

DO  8  J.2,  M  *6 

8BII)>BII)*AII.J|*SIJ)  .7 

RE7URN  ,10 

END  .11 
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SAMPLE  PROBLEM  COMPUTER  OUTPUT 
Rabe's  parameter  =  1.01 

(First  5  lines  are  the  print-out  of  the  input  variables) 


O.«530I535TlO?O*fc*D-O»  0.3040*61  ?*6<>?«92‘I0  00  '0  •  " ri 

0.  T(IS9T1  l*#6e')'»TT60  0?  0.  *6  I  T0*2  9  I  3333  I  ^  ff)  00  0  .  **>636  I  1 6  1 039901  tf-0# 

16  70  ;0  I  1  i  i  0  700  0  3  0  0  0 

0.1000000-00  0. 1000000-00  0. 2000000-00  0.3000000-00  0.2990000  01  1. 0000000-11  -  . 

.0,  -0.  0.1300000  02  0.1000000  01 


T  >  V  0«/nT  OV/OT 

0.  0.30*0*612*6*30  00  -O. • 7*6036380000  00  -0.1298136626*30-01  -0.740*963302910-02 

1.000000000000000  00  0.  0.  0.  -0 

0.  1.000000000000000  00  0.  0.  0 

0.  0.  1.000000000000000  00  0.  -0 

0.  0.  0.  :. 000000000000000  oo  -o 

0.13*3*8300  02  C  3361*770139*0-00  -0.9*38200677*30  00  -0.6105015616370-02  -0.1300066037750-02 
-0.1037197078*96*10  02  -0.000661319*03910  01  0.331031442396930  01  -0.106916302426600  02  -0 

0.263013171619*00  02  0.  1 161*5*'»34  71 760  02  -0.951590100797290  01  0.296411147420660  02  0 

-0.1002810792914*0  02  -0.113129264079120  02  0.095571116454190  01  -0.117062309241000  02  -0 

•0.1*3408020087130  02  -0.614101478204460  01  0.491611141213600  01  -0.162096471074270  02  -0 

0.100471120  02  0.1324436390170-00  -0.929104*121210  00  0.7972371333*00-02  0.1612901049620-02 

-0.1030977043737*0  02  -0.360118300161220  01  0.0111*6930490990  01  -0.100190309429350  02  -0 

0.131936800402090  02  0.123221024411420  02  -0. 1 1609617101' 260  02  0.141011196024240  02  0 

-0.133047832141070  02  -0.110879294013090  02  0.140030007319410  02  -O.11441632201171P  02  -0 

-0.210116128331360  02  -0.811100777416020  01  0.049773309360140  01  -0.21*727723903440  02  -0 

0.462124240  02  0.31U1132797000  00  -C .0138301181390  00  0.1212977301140-01  0.7601346792110-02 

-0.131427734369280  01  -0.924989736341000  00  0.512407063093900  00  -0.172731300164240  01  -0 

0.114209093617460  01  0.326200173109660-01  -0.279003716670710  01  0.2*1302271244400  01  0 

-0.116061998117120  01  -0.219011346316070-00  0.0207139399*4910  00  -0.21*4**01*422*70  01  -0 

-0.201911129743600  01  0.1*8140971011660-00  0.21977100*110140  01  -0.2120124640*3300  01  -0 

0.616766600  02  0.6190146161410  00  -0.7687814403020  00  -0.24201446**130-01  -0.1429116227310-02 

0.1144912120460*0  02  0.10994130910*240  02  -0.91078641*634130  01  0.1204*31**432760  02  -0 

-0.19*414*01632110  02  -0.168136000241140  02  0.170772723240070  02  -0.206112171140010  02  0 

0.1*73706643**110  02  0.1746469974*4110  02  -0.161133444610270  02  0.214216991029230  02  -0 

0.10014072911**00  02  0.0*633007*246670  01  -0.9**3*211413*140  01  0.114*19*77222000  02  -0 

0.771714*10  02  0.3221410491470  00  -0.061*001910790  00  -0.12*9106027*60-01  -0.77*921*311440-02 

0.42*33024111*110  01  0.4241019133710*0-00  -0.284672971166120  01  0.212412141073440  01  -0 

-0.000190097193630  01  -0.143314*41277440  01  0.42127**12434300  01  -0.61344*961311360  01  0 

0.091361173211310  01  0.281209178499260  01  -O.40402001470100O  01  0.343340361414120  01  -0 

0.64***6117*79120  01  0.1*114*722121910  01  -0.297411014946710  01  0.37012*731739070  01  -0 

0.7*3971130  02  0.3040461923230  00  -0.0746*36006660  00  -0. 129*1 117***00-01  -0.74*4930691410-02 

0.119*03022110110-00  -0.1*6227101119470  01  -0.9730143326347*0  00  -0.012290*369127*0  00  -0 
-0.522903374961*90  01  -0.1179*6411*20090-00  0.112047910701000  01  -0.232*31*90721740  01  0 

0.1091**4*6411110  01  -0.240776331941660  01  -0.11293994241*240  01  0.104149747*12*10-01  -0 

J.971662940700120  01  0.116691792762110  01  -0.1140*1349117440  01  0.114**7*913*4260  01  -0 


SU  «*RI*riaN*L  fOU*7ION  CHECKS  -0.237330-11-0.472310-12-0.220460-11-0.173710-11  0.119020-11  0. 
CM*R«C7E*ISTIC  EIPONEN73 

-0.9794309600147*0  00  0.2016*24633*91*0-00  -0.9794309600147*0  00-0.201**24*93*41*0-00 


C 

0.10000712*7137210  01 
.7143700911094440-02 
.1271093461034660-01 
.129*11662644*710-01 
.74*49613*2911010-02 

0.1000071267137210  01 
.71437*0911044040-02 
.1271093461034610-01 
.129*11662644**00-01 
.74*49613*29106*0-02 

0.1000071267137210  01 
.71437*0911094**0-02 
.1271093461034940-01 
.129*1 1**26449060-01 
.74*49613*2914010-02 

0.1000071267137210  01 
.7143700911092910-02 
.1271093461034700-01 

.129*116*2644*640-01 

.74*49613*2912790-02 

0.1000071267137210  01 
.71437*0911094720-02 
.1271093461034740-01 
.129*11**2644*910-01 
.  74*49*13*29126*0-02 

0.10000712*7137210  01 
.71437*0911093200-02 
.1271093461034*60-01 
.129*1 16626449060-01 
. 74*49*13*2911660-02 

0.1000071267137210  01 
.71437*0911091*00-02 
.1271093461034670-01 
.129*116*2644*770-01 
.74*49613*29122*0-02 


176210-11 

T**CE 

0.4109*0799704460-01 


SUO6ES7E0  CH6NGES  IN  INI7UL  VALUES  0.03309194031120-0*  0.7*206031460110-07  0. *21790463*2410-0*  0.344*04442*02*0-07 
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DI/DT 


OV/OI 


0.  0.)0«0«61})1940  00  -0.*T46«%4N1  T<««n  00  -0.  1  2«l  •40^0-01  -0.  74(49 1  T65  )4 1 0-0? 

I .OOOOOOOOOOOOOOD  00  0.  C.  0.  -0. 
0.  1.000000000000000  00  0.  0.  0. 
0.  0.  1. OOOOOOOOOOOOOOD  00  0.  -0. 
0.  0.  0.  1.000000000000000  00  -0. 


O.I)4)4n;7D  02  0. 7  76l4«67(l«60-00  -(' . 94 S B?OSO S?l (0  00  -0. 6l0499^4i^^00-0^  -0.150(0  70144960-0? 


-0. 1(5  719541(91150 
0.?(10149296(6?00 
-0. 1002(07  15  169190 


-0. 

145 

488690)2) 

7  70 

UNOtrioM 

67 

)))02 

IN 

NO 

UNOOIOM 

AT 

12  707 

IN 

NO 

UN0(7iUW 

AT 

12  721 

IN 

NO 

uNO(n  ow 

AT 

))112 

IN 

"0 

UNORFIOM 

AT 

12707 

IN 

NO 

02  -0.8006629992(1770  01 
02  0.  116145714260  190  02 
02  -0.  1  15129490091  no  02 
02  -0.6141024212(4410  01 


0.5510550(2105640  01 
-0.9515929(1114700  01 
0.895571911179710  01 
0.491612864)12(60  01 


-0.1(69164671559(0  02  -0. 

0.29641092796(960  02  0. 
-0.  117062122125250  02  -0, 
-0.162096149926050  02  -0, 


0.10(471090  02  0.1524466262100-00  -0.929194062414000  0.79725472(7140-02  0.1612901101250-02 


-0.1(5897191806760  02 
0.15)956009(595(0  02 


-0.560118(51156670  01 
0.  125221  12249(470  02 


-0.155847072729440  02  -0.110879190776450  02 
-0.2)01156614(4440  02  -0.911189799210950  01 


0.81114756(5198)0  01 
-0.116096459450140  02 
0.148050157)5)250  02 
0.849776172925600  01 


-0.1(0198047202610  02  -0. 

0.  1418)2567969290  02  0. 
-0.1)4415650491900  02  -0. 
-0.21672725501)960  02  -0. 


0.462124260  02  0.510)154282100  00  -0.05(50296(490  00  0.12)2969665410-01  0.7601492008720-02 

-0.151420914026750  01  -0.924946(2)104070  00  0.512442170991800  00  -0.1727464291791)0  01  -0. 

0.  1142  7645501  1470  01  0.525510104(5010-01  -0.2797975)407)620  01  0.2(1489681827(90  01  0. 

-0.  1  160495849445(0  01  -0 .  2 1 9 7 5 (6  112 1 7 200- 00  0.820651679088920  00  -0.21(47552260)140  01  -0. 

0.188191(99)19890-00  0.21976799102951001  -0.2120050)145072001  -0. 


-0.20)90)4(6915140  01 


0.616766660  02  0.619(1)94667)0  00  -1.76(78)998(610  00  -0.2428428107110-0)  -0.1429)056509)0- 


0. 1 1449)404608220 
-0.  19(4  )52(9470190 
0.1875  M002(57210 
0.  10(141001 1 7)6(0 


02  0. 1099414414)9270  02 
02  -0.169156010995060  02 
02  0.  1  74646(60)85240  02 
02  0.886549(40994)20  01 


-0.9)078620916010  01 
0.17077276081)810  02 
-0.1611554)69)2500  02 
-0.9885(290711)100  01 


0.128465)527)1770  02 
-0.2061)248744)840  02 
0.2142172)6)81450  02 
0.1148200705626(0  02 


02 

-0. 

0. 

-0. 

-0. 


0.7717)4810  02  0.522)4289)5410  00  -0.86)9004)27490  00  -0.1269299208010-01  -0.77891825)1460-02 

00  -0.284671552496060  01  0.21242972217)9)0  01  -0. 
01  0.42)27722729)090  01  -0.615444974207720  01  0. 
01  -0.484(261(49667)0  01  0 . 5455(62 194 )5((D  01  -0. 
01  -0.2974296(271 2510  01  0.570126451(6750001  -0. 


0.42(547714187770 
-0.(00 186972594140 
0.(9)556886802  740 
0.648884021(7)970 


01  0.424086745896)80' 

01  -0.1455)212)477970 
01  0.2812061'  1226440 

01  0.  10)47105252950 


0.7(5971150  02  0.50404611)1940  00  -0.8746955817940  00  -0.129(128184050-01  -0.748491765)420-02 

00  -0.6)2)11122)170)0  00  -0. 
01  -0.252848052402190  01  0. 
01  0.10)754509946610-01  -0. 
01  0.1148(5624)8912001  -0. 


0.1)9782421901900-00  -0.1862294^0541280  01 
-0.52290161527)1  in  01  -0.  1  1  796)49590)820-00 
0. )09)R4405581 170  01  -0. 24C7 7(( 7 100 7590  01 
0.5  7166060)8646  711  01  3  .  1  160904U  7  14  1  920  01 


-0.975001924986010 
0. 1)2045708047500 
-0. )12957646)61 150 
-0. 1)40802)2(41 150 


C 

0. )00007)266571410  01 
7 )457)74190)096D-02 
I27)0((l  15018040-01 
129(1281(40529)0-01 
74(491 7651406200-02 

0. 1000071266571410  01 
71457)7419011 190-02 
127)0(8115018040-01 
129(12(184052(70-01 
748491 7651405470-02 


0.100007)266571410  01 
71457174190)2)40-02 
127)0881 1501 79)0-01 
129(12(184052960-01 
74(491 7651404290-02 

0.100007)266571410  01 
71457)74190)2790-02 
127)08(1  15317950-01 
129(12(184052960-01 
748491 755)404280-02 

0. )00007)2665714in  01 
7)4»7)74190))89D-02 
127)08(115017(80-01 
129(12(18405)050-01 
748491765)40)))D-02 

0. )00007)266571410  01 
7)457)74190)4120-02 
127)08(115018090-01 
129(12(18405)190-01 
74(491765)405110-02 

0.100007)266571410  01 
7)457)74190)1180-02 
127)08(115017920-01 
129(12(1(4052960-01 
748491765)40)920-02 


S  *  V6RI6t|0NAl  fOOAIION  CHfCKS  -0.191490-11  0.1065(0-12-0.266100-11-0.116620-11  0.588420-12  0.155790-H 
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CHAHACTfHISTIC  EIPONfNTS  I«*Cf 

•0.9r94906««e6jl*n  oo  0.;0l6fl)<)a6)7(in<i0-00  -0.9T9450696a62 1*0  00-0.20l6839a6JT8090-00  0.91098  r062T3T220-0l 

SUCCE STEP  CHANCES  IN  INIII81  » *tUt S-0. 289 T012 1  8 T91  20-0 j- 0. 1 6681 029 1 T900n-0»  0 .  I  6  I T0899626 1  80-0 ) -0 . 28 1 T22 T929 1  6  TO- 0  J 
SUGGESTED  CHANGES  IN  INITIAL  W AlUl S -0. i8 10 1026 T 10 1 TO-0 T-0 . 2  I  969992968 T20-0f  0 . 2  1  96 1 9 1 68 T S800-0 ? -0 . 1 T 16809 19 TS220- 0 T 
NEM  GUESSES 

0.9090960990908060  00  -0. 8T9689601 T6ia89li  00  0.  86  I T09 121 9 198 190  00  0.996961  1781221810-00 


81ANF  arSTAICTEO  tH8EI  800T  MUBIEM 
SUN  JUPHIA 

8N  •  0.99187919710708680-01 

JACnai  CUNSTANI  •  O. 100007126697|9|00  Ol 

PEAIUO  •  0.78997119868997760  02 

«  VP  SUB  A 

initial  values  0.90909611119181210  00  -0.87968998179199660  00  0.86170929999191790  00 

FINAL  VALUIS  0.90909611119996870  00  -0.87968998179198210  00  0.86170929999109660  00 

INITIAL  0«/0I  ANO  Ov/OT  VALUES  HfHI  -0.12981281890929270-01  -0.79899176919061980-02 

ONE  OF  THE  CONJUOAIE  CONPlTA  C MA8 AC T f M I S T I C  ROOTS  IS  -0.97999069686211910  00  0.20168198617809010-00 

(Ml  TRACE  mas  0.91098706279721990-01 


P  SUB  V 

0.99696121999092610-00 

0.99696121999109610-00 


